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ABSTRACT 

We develop a numerical code to calculate the neutrino transfer with multi- 
energy and multi-angle in three dimensions (3D) for the study of core-collapse 
supernovae. The numerical code solves the Boltzmann equations for neutrino 
distributions by the discrete-ordinate (S„) method with a fully implicit differenc- 
ing for time advance. The Boltzmann equations are formulated in the inertial 
frame with collision terms being evaluated to the zeroth order oi v/c. A basic 
set of neutrino reactions for three neutrino species is implemented together with 
a realistic equation of state of dense matter. The pair process is included ap- 
proximately in order to keep the system linear. We present numerical results 
for a set of test problems to demonstrate the ability of the code. The numerical 
treatments of advection and collision terms are validated first in the diffusion and 
free streaming limits. Then we compute steady neutrino distributions for a back- 
ground extracted from a spherically symmetric, general relativistic simulation of 
I5M0 star and compare them with the results in the latter computation. We 
also demonstrate multi-D capabilities of the 3D code solving neutrino transfers 
for artificially deformed supernova cores in 2D and 3D. Formal solutions along 
neutrino paths are utilized as exact solutions. We plan to apply this code to the 
3D neutrino-radiation hydrodynamics simulations of supernovae. This is the first 
article in a series of reports on the development. 

Subject headings: methods: numerical — neutrinos — radiative transfer — stars: 
massive — stars: neutron — supernovae: general 



- 3 - 



Introduction 



Elucidating the explosion mechanism of core-collapse supernovae is a grand challenge 
in astrophysics, nuclear and particle physics as well as computing science. It requires 
the numerical computations of hydrodynamics, neutrino transfer and electromagnetism 
that incorporate detailed microphysics at extreme conditions. Delicate interplays of these 
multiple elements in various physics have to be described quantitatively in numerical 
simulations in order to determine unambiguously the outcome of collapse and bounce of 
central cores and to obtain the magnificent optical display, neutrino burst, and ejected 
material as well as remnants that fit observations appropriately. Despite the extensive 
studies with ever increasing theoretical achievement and computational resour ces for more 
than four decades, the critical e 



Kotake et al. 



2006 



Janka et al. 



emen t of supernova explosion has been elusive (iBethe 



1990 



20071 ). 



One of the key issues is the neutrino transfer. Interactions of neutrinos with material 
play important roles in a couple of ways. The lepton fraction stored in the bouncing core is 
determined by the neutronization and neutrino trapping; the core bounce launches a shock 
wave, which stalls on the way out owing to energy losses through the dissociation of nuclei 
and the neutrino cooling; the fate of the stagnated shock wave is eventually determined by 
the flux of neutrinos that are emitted copiously by a proto-neutron star and interact with 
the material below the shock wave. Neutrino signals detected at terrestrial detectors are 
an important probe into what is happening in th e deep interior of t hickly veiled supernova 



cores, which was indeed vindicated for SN1987A f lHirata et al 



19871). 



Among viab 



teethe fc Wilson 



e scen arios for supernova explosions, the neutrino heating mechanism 



19851 ) is currently the most promising. The heating of material below 
the stalled shock wave by the absorption of neutrinos emitted from deeper inside the core 
is supposed to revive the stagnated shock wave. Since neutrinos carry away most of the 



-4 - 



liberated energy of ~10^^ erg, tapping a hundredth the energy transferred by neutrinos is 
sufficient to obtain the observed explosion energy of ~10^^ erg. Even if the neutrino heating 
is not the main cause of explosions, it will be still critical to accurately evaluate the energy 
exchanges between neutrinos and matter. For this purpose, one has to solve the neutrino 
transfer equations and obtain fluxes, energy spectra, emissions and absorptions. 

The numerical treatment of neutrino transfer is a longstanding problem though. 
One of the difficulties is the wide range of opacity in the supernova core. Deep inside, 
interactions are so frequent that thermal and chemical equilibrium is established in a 
quasi-static manner and neutrinos spread in a diffusive way. As the density and temperature 
become lower in the outer region, the neutrino interactions become less frequent and the 
statistical equilibrium is not achieved. Further outside, neutrinos propagate freely through 
transparent material. It is the intermediate regime between diffusion and free-streaming 
that the neutrino heating mechanism operates. It is hence mandatory in the computation 
of neutrino transfer to treat all these regimes. 

It is computationally demanding, however, as one often sees in the applications of 
radiation transfer to other subjects. Even in spherical symmetry, the neutrino transfer is 
a three-dimensional problem (or four-dimensional if one includes time) since the neutrino 
distribution is a function of neutrino energy, angle of neutrino momentum with respect 
to the radial direction as well as radial position. Since the neutrino interactions are 
strongly energy- dependent, the multi-energy group treatment is indispensable. Without 
any symmetry, the neutrino transfer becomes a six dimensional problem: the neutrino 
distribution is then as a function of three spatial coordinates and three neutrino momentum 
(energy and two angles that specify the direction of neutrino momentum are commonly 
used as independent variables). With limited computational resources, this fact makes the 
numerical computation of neutrino transfer such a formidable task and various levels of 
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approximations have been employed (See ^for details) so far. 



The numerical treatment of neutrino transfer under spherical symmetry has been 
sophisticated over the years - so much so that the Boltzmann equation is solved without 
any approxim ation. Unfortunately, it has been consisteiitly demonstrated with scrutini zed 



microphysics (IRampp fc Janka 



Thompson et al. 



2003 



2000; 



Liebendorfer et al. 



Sumiyoshi et al 



2001 



Mezzacappa et al. 



2001 



20051 ) that no explosion occurs under spherical 



symmetry. On the other hand, recent multi-dimensional simulations have revealed the 
critical role of hydrodynamical i nstabilities with/without possible neutrino assistance for 



successful supernova ex plosions (IBurrows et al. 



2010 



Suwa et al. 



2006a : 



Marek fc Janka 



2009 



Bruenn et al. 



20101 ). Core-collapse simulations in axial symmetry have adopted 



different approximations and only recently the Boltzmann equations wer e directly solvec 



by the discrete-ord inate method for a handful of post-bounce evolutions (lOtt et al 



Brandt et al. 



2008 



20111 ). In three spatial dimensions (3D) without any symmetry, most of 



the simulations done so far focused on the hydrodynamical instabilities with very simple 
treatments of neutrinos. Very recently the first 3D simulations with the neutrino transport 
being treated with th e isotropic diffusion source approximation have been performed 



(ITakiwaki et al. 



20111 ). Unfortunately, their results are inconclusive owing to limited 
numerical resolutions. 3D computations of neutrino transfer have just begun and further 
development is highly awaited. 

Exploring the explosion mechanism in three dimensions is indeed the current focus 
in the supernova society. This is mainly because marginal explosions have been observed 
in many 2D simulations so far and it is expected that 3D will give a boost: increased 
degrees of freedom may give fluid elements more time to hover arou nd the heating region . 



It was indeed demonstrated by a systematic numerical experiment (iNordhaus et al 



20101) 



that the critical neutrino luminosity is lowered and explosions occur more easily as the 
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spatial dimension increases (See also iHanke et al.ll201l[ ). This needs confirmation by more 
realistic simulations with detailed neutrino transfer. It is interesting to assess quantitatively 



whether the observed rotation and/or magnetic fields of pu 



hydro dynamica . 



Rantsiou et al. 



sars can be induced by 3D 



inst abilities as suggested by some authors ( iBlondin fc Mezzacappal 12007 
20111). Moreover, gamma ray bursts may be produced by anisotropic 

2010al Jbl). 3D neutrino transfer in these contexts has not 



neutrino radiations (Harikae et al. 



been studied yet and will be the target in our investigations to come. 

Under these circumstances, we have developed a numerical code to treat the time- 
dependent neutrino transfer in 3D space with no symmetry. We adopt a discrete-ordinate 
method (S^) and deploy multi-angle, multi-energy groups to directly solve the Boltzmann 
equation for the neutrino distribution function in six dimensional phase space. The fully 
implicit finite differencing in time is employed. The neutrino reactions and EOS of dense 
matter that are suitable for supernova simulations have been implemented. 

As far as the authors know, this is the first study on the 3D neutrino transfer for 
core-collapse supernovae. The purpose of this article is to report the first results on the 
performance of the new 3D code. As a first step, we study the neutrino transfer in static 
backgrounds in this paper. Starting with some basic tests in idealized settings, we then 
examine the code performance for exemplary core profiles both before and after bounce 
with the appropriate microphysics inputs. By these tests and applications, we demonstrate 
that the 3D simulation of neutrino transfer is now feasible in all the opacity regimes. It is 
true that the problem size is very large and the resolution is still severely limited by the 
current computing resources. We discuss the computing power required for the full 3D 
core-collapse simulations based on the current numerical experiments. 

We arrange the article as follows. We briefly describe the recent developments of 
neutrino transfer computations for supernovae in [|2l In [|3], we give the basic equations of 
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neutrino transfer together with the microphysics included in the current version of the code 
and explain their numerical implementation. We report in ^the results of the basic tests 
to validate the code: the neutrino transfer in the diffusion and free-streaming limits are 
tested in §4.1( the implementations of neutrino reactions as collision terms are examined in 
§4.21 After these basic tests, we proceed to some applications of the code to more realistic 
background models adopted from supernova simulations in §3 we first compute spherically 
symmetric neutrino transfers in §5. II and compare them with the results published in other 
papers. We then demonstrate some basic features of 2D and 3D neutrino transfers using 
artificially deformed core profile in § §5.21 and 15.31 The summary is given in §H] with some 
discussions on further extensions of the current code. 



2. Developments on Neutrino transfer 



We describe briefly the recent developments of neutrino transfer in core-collapse 
supernovae in order to grasp the status of prescriptions to handle the neutrino transfer in 



the previous studies. Further re: 
we refer the rey iew articles by 



erences to overview the historical and modern deve^ 



Qtt et al. 



Suzukil (119941 ): 



Kotake et al. 



(|2Q06|); 



(120091 ) on core-collapse supernovae as well as the books by 



Janka et al. 


( 


2007 


); 


Peraiah 


(2002 


); 


Castor 



(120041 ) on the radiation transfer. 



Under spherical symmetry, the development of neutrino transfer has reached a level 
of sophistication by direct solutions of equations, thanks to enough computing resources 
recently. After longs tanding effor ts on the treatment of neutrino transfer by invoking 



approximations (See 
transfer with mu 



Suzuki 



( Rampp &: Janka 



19941 . for example), the numerical solution of the neutrino 
ti-energy and multi-angle has become possible under spherical symmetry 



2000 



(ILiebendorfer et al 



2001 



lompson et al 



2003 



Liebendorfer et al. 



Mezzacappa et al 



2004 



Sumiyoshi et al 



20011) in general relativity 



2005|). Through the 
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developments of the frameworks, the method of numerical solutions both in inertial and 
co-moving fram es have been estab lished to handle the neutrino transfer with the Lorentz 



transformation (IBuras et al 



20061 ). The comparison of the methods to solve the neutrino 



transfer by direct solutions of the Boltzmann equations and by iterations of the moment 



equations with a variable Eddington 



of neutrino-radiation hydrodynamics (iLiebendorfer et al. 



actor has been also mad e to test the approaches 



20051 ). The microphysics of the 



equation of state and the neutrino reactions have been implemented to examine their 
influence in detail by having the advance of nucle ar physics in superno vae. The most 



updated set of neutrino reactions in dense matter ( 



Burrows et al. 



capture rates of nuclei (ILanganke fc Martmez-Pinedc 



impa ct on the explosion mechanism ( 



Langanke et al. 



2003 



2006bl ) as well as electron 



2003 ) have been adopted to test t 



Hix et al. 



2003 



ne 



Buras et al. 



20061 ). The role of the equation of state based on the unsta ble nuclei has been exa mined 



20051 ). Even 



through the delicate connections with the neutrino transfer (iSumiyoshi et al. 
with the exact treatment of the neutrino transfer and the updated sets of microphysics, the 
numerical studies of the gravitational collapse of massive stars have shown that there is no 
successful case of supernova explosion under spherical symmetry for sets of stellar models. 

In the developments of the numerical methods of neutrino transfer under spherical 
symmetry, there has been progress on clarifying the role of neutrino transfer in core-collapse 
supernovae. The proper treatment of neutrino transfer is crucial to determine the amount 
of neutrino trapping in the collapsing core, the energy and flux of neutrinos emitted 
from the cooling region, the heating behind the stalled shock wave and the prediction 
of the spectrum of supernova bursts. The progress from the early methods such as the 
light-bulb approximation, the leakage scheme and the flux limited diffusion method to 
the exact treatment of neutrino transfer has clarified the influence of approximations to 
those quantities. Among others, the neutrino heating is sensitive to the neu trino transfer. 



especially, to the flux factor of neutrino distributions f lJanka &: Mueller 



19961). The accurate 
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evaluation by the neutrino transfer is necessary since the flux limited diffusion method ma y 



overestimate the flux factor, resulting underestimation of the heating ( lYamada et al. 



19991 ) ■ 



It is to be noted that the proper treatment of neutrino transfer enables ones to predict 
the energy spectrum of the supernova bursts for terr estrial observations (IThompson et al. 



2003 



Sumiyoshi et al. 



2005 



2007 



Fischer et al 



20091). The energy spectra predicted by the 



neutrino-radiation hydrodynamics are used to evaluate the event numbers for supernova 
explosions and black hole formations near the Galaxy in future by t aking into account 



the neutrino oscil 



Ando et al. 



2005 



ations and the speci f ication of neutrino detect ors (ITotani et al 



Nakazato et al. 



2010 



Keehn &: Lunardini 



19981: 



201 



In two dimensions, the approximate treatments of the neutrino transfer have 
been adopted in most of the recent studies by state-of-the-art calculations. This is 
because the exact treatment of neutrino transfer becomes much harder than the case of 
spherical symmetry due to the increase of dimension of phase space from three to five. 
There are several categories in the approximate methods, depending on the degree of 
accuracy regarding transfer and dimensional assumption. Putting the emphasis on the 
two dimensional behavior, the flux limited diffusion method has been adopted for the 



neutrino-radiation hydrodynamics to rey eal the explosion mechanism (ILivne et al. 



2004 



Walder et al. 



2005 



Burrows et al. 



2006aJ). The diffusion equations of neutrino energy 
and flux distributions with multi-energy groups are solved with the flux-limiter to handle 
the transition to the free-streaming limit. Although this approach is advantageous to 
describe the lateral transport of neutrinos in supernova cores, the intermediate regime from 
diffusion to transparent regimes are handled by the prescribed flux-limiter. Attempting 
to describe better the transition of neutrino transfer in the intermediate regime, the 
method of "ray- by-ray plus" approximation has been ad opted for the neutrino-radiation 



hydrodynamics ( iBuras et al 



2006 



Marek fc Janka 



neutrino transfer for spherical symmetry (IRampp fc Janka 



20091) utilizin g the developed code of 



2OO2I ). The equation of neutrino 
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transfer along each radial ray is solved independently for many directions to cover the whole 
region of the supernova core. By using the exact treatment of the ID neutrino transfer for 
multi-energy groups, it is advantageous to describe the whole regime of neutrino transfer. 
By assuming the spherical symmetry on the transport along each ray, the transport of 
neutrinos for lateral directions is neglected except for the advection with material and a 
partial contribution of pressure. This prescription along independent rays overestimates 
the angular dependence of neutrino quantities and enhances the neutrino fluxes along the 
radial directions. Some of the recent studies adopt the mixed approach of th e ray-by-ray 



plus 



approximation together with the approximation of flux li mited diffusion ( 



2006f) or the isot ropic diffusion source approximation (IDSA) (ILiebendorfer et al. 



Bruenn et al. 



2009 



Suwa et al. 



20101 ). The ray-tracing met hod has been utilized for the analysis of anisotropic 



neutrino radiations ( iKotake et al. 



20091). 



Recently, the 2D numerical simulations by the multi-angle, multi-energy group 
neutrino-radiation hydrodynamics have been done to study the postbounce phase of 
core-collapse supernovae. The neutrino transfer is solved by the discrete-ord inate (S„) 



20041 ). In order 



method to describe the whole regime of neutrino transfer (iLivne et al. 
to save computational load and to have large time steps, the flux limited diffusion 
approximation is adopted fo r the central pa rt of supernova core for the study on the time 



evolution for a long period (lOtt et al 



20081 ). The basic behavior of neutrino quantities in 



2D has been reported through the examination of angle dependence, moments of energy and 
angle of neutrino distributions in realistic profiles of supernova dynamics. The advantage 
to describe the neutrino transfer in 2D has been demonstrated through comparisons with 
the counterpart by the flux limited diffusion approximation. The multi-angle treatment of 
neutrino transfer has revealed the substantial effects such as asymmetries in neutrino fluxes 
between pole and equator, enhance ments of the neutrino heating throu gh the integral of 



neutrino sources over many angles (lOtt et al 



2008 



Brandt et al. 



201 If ). With those levels 
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of approximations and sophistications described above, the mechanism of core-collapse 
supernovae has been explored by neutrino radiation hydrodynamics in 2D to find a handful 
of successful models with new mechanisms, which are still under debate. 

In three dimensions, the trea tment of neutrino transfer is still in its infancy. In most o: 



the numerical simulations in 3D ( iBlondin et al. 



2003 



Ohnishi et al. 



2006 



Blondin fc Shaw 



2007 



Iwakami et al. 



20081 ) . the neutrino transfer has been ne glected or simp 



assuming the gjiven ^ rate of neutrino cooling and/or heating ( iNordhaus et al 



ified 



by 



2010 



Rantsiou et al.ll2011i) . The central part of supernova core was often omitted to avoid 
the treatment of neutrino trapping and emission in the early stage of researches. These 
simplifications in the numerical studies are made to explore the frontier of hydro dynamical 
instabilities and to seek the favorable conditions of explosion in 3D beyond the assumption 
of axial symmetry. Moreover, the computational cost of neutrino transfer in 3D has 
been formidable to handle the time evolution of neutrino distributions in six dimensions. 
The 3D numerical simulations by the ray-by-r ay plus approximat ions have been reported 



recently by adopting the fiux li mited diffusion 



source approximation (IDSA) (ITakiwaki et al. 



Bruenn et al. 



20101 ) or the isotropic diffusion 



201ll ) for radial transport. However, the 



direct solution of the neutrino transfer in 3D has not been implemented in the numerical 
simulations so far. 

There have been long-term efforts on the numerical method to solve directly the 
neutrino transfer in multi-dimensions besides the approximative approaches mentioned 
above (See also. 



Swesty fc Myra 



relativity have been studied by 



2009). General forms o 



Cardall fc Mezzacappal ( l2003l ) 



the n eutrino transfer in g eneral 



Cardall et al. 



(120051). A 



time-dependent Boltzmann transport scheme in multi-energy and multi-angle has been 
recently developed for neutrino-radiation hydrodynamics in one and two dimensions 



(ILivne et al. 



20041 ). The Boltzmann equation in 2D axisymmetric geometry is discretized 



- 12 - 



in conservative form using the discrete-ordinates (S„) method by dropping the velocity 
dependent terms. The 2D transport method incorporated to the neutrino radiation 
hydrodynamics is apphed to a time-dependent 2D test of a post-bounce supernova core. 
The neutrino radiation hydrodynamics code (VULCAN-2D) has been utihzed for the 
2D supernova simulations with a va riant of the flux limited method as mentioned above 



dOtt et al. 



2008 



Brandt et al.li2011i) . More recently, a new algorithm to solve the neutrino 



transfer in two dimensions has been developed to conform the Lorentz transformation in 



the transport equation (IHubeny &: Burrows 



20071 ). They derive the formulation using the 



mixed-frame approach by evaluating the collision term in the comoving frame with a Taylor 
expansion regarding Lorentz shifts. The new formalism has been applied to one-dimensional 
tests of stationary solutions and proto-neutron star cooling. 

Our study here is to establish the numerical solver of the Boltzmann equation for 
neutrinos in three dimensions, for the first time, beyond the previous developments in two 
dimensions. We develop a numerical code to solve the Boltzmann equation for multi-energy 
and multi-angle group in 3D spatial coordinates. We take an approach to solve the 
Boltzmann equation in the inertial frame, on which we report below, as a basis for our 
developments. We extend the formulation and its numerical implementation by evaluating 
the collision term according to the Lorentz transformation as a next step, which will be 
reported separately elsewhere. 



3. Formulations 

3.1. Boltzmann Equation 

In our numerical code for the neutrino transfer, we solve the Boltzmann equation for 
the neutrino distribution by a discrete ordinate (S„) method. Our starting point is the 
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Boltzmann equation, 



1^ + ^ 

c dt ds 



15/ 

c 5t 



(1) 



collision 

for the neutrino distribution function, f{r,t;e,n), at position, r, and time, t, along path 
length, s. The right hand side is the collision term, which expresses the time rate of change 
due to the neutrino reactions such as emissions, absorptions and scatterings. We prepare 
the neutrino distributions, 

r(r,^,0,t;£-,n-), (2) 

where e*" and n*" are the neutrino energy and the unit vector of neutrino momentum, 
respectively, in the inertial frame. We adopt spatial variables, r, 6, (p, in the spherical 
coordinate system. The unit vector of neutrino momentum is defined with respect to the 
radial direction along the coordinate r as in Fig. [1] We adopt the neutrino angles, 6^, (f)^, 
and the neutrino energy, e*", to designate the neutrino momentum. 

We take an approach in the inertial (laboratory) frame to write down the equation of 



neutrino transfer and to handle the neutrino quanti ties. The way of so 



transfer differs very much depending on the frame (iMihalas &: Mihalas 



utions of neutrino 



19991). The two 



major ways in the comoving and inertial frames have both easiness and difficulty in the 
procedures of solution. On the one hand, the form of left hand side of Eq. ([1]) is simple 
in the inertial frame, while the derivative terms i n the left hand sid e are complicated with 



20061 ). On the other hand. 



velocity dependent terms in the comoving frame (iBuras et al. 
the collision term can be calculated easily in the comoving frame, where the neutrino 
reactions occur in the moving fluid. The collision term in the inertial frame requires 
tedious procedures through the Lorentz transformation of reaction rates from the comoving 
frame in principle. In our strategy, we take the simplicity of the transfer equation in three 
dimensions and will make numerical efforts to handle the collision term in a next step of 
the development. 
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Fixing the framework in the inertial frame, the Boltzmann equation, Eq. ([T]), in the 
spherical coordinate system is expressed as 



1 dr ^ . df 
+ cos 



c dt 
sin^ df 



^ sin 9y cos (9/*" ^ sin 9,j sin 



(9cos 6*^ 



(9r ' r 86 

sin 6',^ sin cos 6 df 
r sin 9 deb 



rsin 6' 



90 



with the definition of the neutrino direction angles (fomraning 



(3) 



•miiision 



19731 ). We remark that 



there is neither velocity-dependent term nor energy derivative in the equation in the inertial 
frame, being different from that in the comoving frame. Choosing the angle variable 
IXy = cos 6iy instead of 6^, the equation can be written by 

1 dp"" 9/'" fJ'l cos (pu df" fil sin (p^ df" 



c dt ~^ dr 
^l-p^9f« 
r d^iy 



de 



1 — /i^ sin cos 6 df 
r sin 9 dd>. 



rsin 6 

16f''' 



c St 



(4) 



collision 



For the numerical calculation, we rewrite the equation in the conservative form as. 



Idf^ _^ ^ A(^2y.^n^ ^ ' 

c dt dr rsin 9 



cos 



'^^sin er) + ^^ 



sm 



df 



de 



■\/l — /i^ cos 9 d 



sm 



rsin 9 

'16f" 

c~5r 



(5) 



- collision 



r d^ifj ^ r sin 9 

We adopt this equation as the basis for our numerical code. We remark that the neutrino 
distribution function is a function of time and six variables in the phase space as written by, 



r{r,9A,t■.^^uAu.e'"). 



(6) 



In the above expressions, the angle variables, iiy and (pi, are those measured in the inertial 
frame. 



3.2. Neutrino Reactions 



We implement the rate of neutrino reactions with the composition of dense matter 
as contributions to the collision term. We take here several simplifications to make the 
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neutrino transfer in 3D feasible. 

As the first step of 3D calculations, we treat mainly the case of static background of 
material or the case where the motion is very slow so that f/c is very small. In the current 
study, we evaluate the collision term of the Boltzmann equation to the zeroth order of f/c 
by neglecting the terms due to the Lorentz transformation. For dynamical situations in 
general, this drastic approximation will be studied carefully by evaluating the effects from 
the Lorentz transformation in future. We plan to implement such effects in all orders oi v/c 
in our formulation by taking into account the energy shift by the Doppler effects and the 
angle shifts by the aberration in the collision term. 

In addition, we limit ourselves within a set of neutrino reactions to make the solution 
of Boltzmann equation possible in the current supercomputing facilities. In order to avoid 
the energy coupling in the collision term, we do no t take into account en ergy- changing 
scatterings such as the neutrino-electron scattering (IBurrows et al.ll2006al ). This makes the 
size of the block matrix due to the collision term smaller and the whole matrix tractable in 
the system of equations. As a further approach, we linearize the collision term for the pair 
process to avoid the non-linearity in equations and to guarantee the convergence. 

In future, having enough supercomputing resources, we will be able to include the 
energy-changing reactions by enlarging the size of block matrices. We also may be 
able to solve the full reactions by the Newton iteration, which requires the complicated 
matrix elements by de rivatives, as have been accomplished in the spherical calculations 



( ISumiyoshi et al 



20051) 



In the numerical study under the assumptions above, we implement the collision term 
in the following way. We utilize directly the neutrino distribution function in the inertial 
frame to evaluate the collision term. We use the energy and angle variables in the inertial 
frame in the calculation of the collision term by dropping the shifts. We drop the superscript 
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in for the inertial frame in the following expressions. For the emission and absorption of 
neutrinos, the collision term for the energy, s, and the angles, fii, and (p^y, is expressed as. 



c St 



(7) 



emis—abs 

Hereafter we suppress the spatial variables and use Q to denote the two angle variables for 
the compactness of equations. The emission rate is related with the absorption rate through 
the detailed balance as, 

R^rnisie, n) = R,,,{e, n)e-^^'-'^-^\ (8) 

where (3 = l/ksT is the inverse of temperature and = iip + fi^ — A*n is the chemical 
potential for neutrinos. The collision term for the scattering is expressed by, 

scat . 

de'e''' 



iSf 
c St 



+ 



J dn' R,,,t{e,n;e',n')f{e,n)[l- f{e',n')] 



(9) 



where fl' denotes the angle variables after/before the scattering. The energy integration 
can be done by assuming the iso-energetic scattering. The expression can be reduced as 

^2 



is_i 

c St 



scat 



(2vr) 



dn' R,,at{^;n')[f{e,n)-f{e,n')], 



(10) 



with the relation, Rscati^'',^) — Rscati^',^')- The collision term for the pair-process is 
expressed by. 



-isf 




r de'e'^ f 


c St 


pair -J 


1 {2nf j 



dQ' R. 



■pair—anni 



+ 



de'e 



dQ' Rpair-emis{e,n;e',n')[i - f{e,nm - f{e',n% 



(11) 



where f{e', fl') denotes the distribution of anti- neutrinos. From the detailed balance, the 
following relation holds; 



^pair—anniy'^ J ii, c , j — iXpair—emisy'^j Ji, t , ye 



(12) 
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We linearize the collision term, Eq. ffTTj) . by assuming that the distribution for anti-neutrinos 
is given by that in the previous time-step or the equilibrium distribution. This is a good 
approximation since the pair-process is dominant only in high temperature regions, where 
neutrinos are in thermal equilibrium. We adopt the approach with the distribution in the 
previous time-step in all of the numerical calculations with pair processes in the current 
study. We utilize further the angle average of the distribution when we take the isotropic 
emission rate as we will state. We have also tested that the approach with the equilibrium 
distribution determined by the local temperature and chemical potential works equally well. 



As for the reactio n rates, we take main 



with some extensions ( ISumiyoshi et al. 



20051 ) 



y from the conventional set by iBruennI (119851 ) 
. We implement the neutrino reactions. 



for the absorption/emission, 



e~ + p i — > z/e + n [ecp] , 
+ n i — y i^e + P [aecp] , 
e" + A i — > z/g + A' [eca] 

v + Ni — >v + N [nsc], 
u -\- A i — > u + A [esc] , 



(13) 
(14) 
(15) 

(16) 
(17) 



for the iso-energetic scattering. We do not take into account the neutrino-electron scattering. 
It is well know n that the influence of this reaction is minor although it contributes to the 



thermalization (IBurrows et al. 



2006al ). As for the pair-process, we take the electron-positron 



process and the nucleon-nucleon bremsstrahlung as follows, 

+ i — > Ui + Vi [pap] , 
N + Ni — > N + N + + [nhr]. 



(18) 
(19) 



For these pair processes, we take the isotropic emission rate as an approximation, which 
avoids complexity but describes the essential roles. We remark that the set of the reaction 
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rates adopted in the current study is the minimum, which describes sufficiently the major 
role of neutrino reactions in the supernova mechanism. Further implementation of other 
neutrino reactions and more sophisticated descri ption of reaction rates in the modern 



version (IBuras et al. 



2006 



Burrows et al. 



2006bl ) will be done once we have enough 



computing resources in future. 



3.3. Equation of State 



We utilize the physical equation of state (EOS) of dense matter to evaluate the rates 
of neutrino reactions. It is necessary to have the composition of dense matter and the 
related thermodynamical quantities such as the chemical potentials and the effective mass 
of nucleon. We implement the subroutine for the evaluation of quantities from the data 



table of EO S as used in the other simulations of core-collapse s upernovae (iSumivos 



2005 



20071 ). We adopt the table of the Shen equation of state (IShen et al. 



1998a 



li et al. 



20111 ) 



in the current study. Other sets of EOS can be used by simply replacing the data table. 



3.4. Numerical Scheme 



We describe the numerical scheme employed in the numerical code for the neutrino 
tra nsfer in three dimensions. Th e method of the discretization is based on the approach 



Mezzacappa fc BruennI ( 



Swesty fc Myra 



(120091) 



1993); 



Castor 



Stone et al 



(120041 ). We refer also the references by 



( I1992I ) for the other methods of discretization of 



neutrino transfer and radiation transfer. 

We define the neutrino distributions at the cell centers and evaluate the advection 
at the cell interfaces and the collision terms at the cell centers. We describe the neutrino 
distributions in the space coordinate with radial A^^-, polar Nq- and azimuthal A^^^-grid 



- 19 - 



points and in the neutrino momentum space with energy A^g-grid points and angle Nq^- and 
A^0^-grid points. We explain the detailed relations to define the numerical grid in §A.2[ 

We discretize the Boltzmann equation, Eq. ([5]), for the neutrino distribution, in a 
finite differenced form on the grid points. Here we assign the integer indices, n and n + 1, 
for the time steps and, i, for the grid position. We adopt the implicit differencing in time to 
ensure the numerical stability for stiff equations and to have long time steps for supernova 
simulations. We solve the equation for /""'"^ by evaluating the advection and collision terms 
at the time step + 1 in the following form, 

71+1 



c At 



+ 



n+1 



1 — COS (py d 



36 



+ 



1 — /i^ sin df 



rsin 9 



rsin 9 

n+l 

+ 



;sin 9f) 



1 d 



n+l 



fj^l COS 9 d 



sin 9 



(sin (pyf) 



- n+l 






_c 6t_ 



n+l 

collision 



(20) 



where we schematically express the advection terms for the cell containing fj^~^^. We 
evaluate the advection at the cell interface by the upwind and central differencing for 
free-streaming and diffusive limits, respectively. The two differencing methods are smoothly 
connected by a weighting factor in the intermediate regime between the free-streaming 
and diffusive limits. We describe the numerical scheme for the evaluation of the advection 
terms in §A.3[ We express the collision terms by the summation of the integrand using the 
neutrino distributions at the cell centers. 



3.5. Solution of Linear Equation 

We arrange the discretized neutrino distribution as a vector for the system of linear 
equations. The length of the vector is N^ector = NrNsN^NeNs^N^^. We advance time-step 
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by the implicit method for a time step, At, by the relation, 

for the vector of the neutrino distribution, Z""*"^, at the time step, n + 1. By linearizing the 
collision term as described above, we rearrange the Boltzmann equation as a set of linear 
equations, Af^^^ = b. The matrix. A, contains the terms from time-advance, advection 
and absorption terms. The source vector contains the old vector and emission terms. 
This large sparse matrix {N^ector x N^ector) contains block diagonal matrices of the size, 
Ng^N^^ X Ng^N^^, which comes mainly from the angle- coupling due to scatterings, together 
with the lines of none- zero elements due to the advection in the three directions. We remark 
that the size of the block matrices becomes huge as N^^Ng^N^^ x N^Ng^N^^, if we take the 
energy changing reactions and the Lorentz energy shifts fully into account. 



(iSaad 



We sol ve this system of equations by the matrix solver using the iterative method 



20031). We use the B i-CGSTAB method by utilizing a program in the Templates 



library fiBarrett et al. 



19941 ) with the point- Jacobi method as a pre-conditioner. We set 
the allowable convergence measure to be 10~® and get the convergence typically within 20 
iterations for the current numerical studies. 

We solve the evolution of neutrino distributions for multi-species (Nf). For the basic 
tests, we treat the one species of neutrino {Nf = 1). For the applications of supernova 
cores, we treat the three species of neutrinos, z/g, and z/^ {Nf = 3). We treat z/^ as a 
representative of the group for four species z/^, u^, Ur and Ur- 



4. Basic Tests 



We performed the series of basic tests on the advection and collision terms, in order 
to validate the numerical code to solve the neutrino transfer in three dimensions. We 
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report here the two tests on the advection in the diffusion and free-streaming hmits, where 
we can compare with the analytic behavior. These tests are performed to vahdate the 
advection part of the Bohzmann solver. In addition, we report the tests on the stationary 
solution, the time evolution toward the equilibrium and the comparison with the spherical 
calculations. These tests are used to check the collision term as a source term and to 



examine the neutrino reactions with dense matter in supernova cores. T 



le method of basic 



tests in the current s tudy are based on the standard tests described in (IStone et al 



Swesty fc Myra 



1992 



20091). 



4.1. Advection Term 

4.1.1. Diffusion Limit 

In order to show that the new code can properly handle diffusions of neutrinos in 
opaque material, we compute the diffusion of a Gaussian packet in a uniform background. 
Taking a scattering as a sole contribution to the collision term, we assume that it is 
isotropic and isoenergetic and its rate is independent of incident neutrino energy. We keep 
multiple energy bins to confirm that the neutrino energy spectrum is unchanged during the 
simulation. 

Analytic solutions are available for this test. The diffusion of a Gaussian packet with 
its central position located at rg and its width, d^, being (ADtoY^'^ initially is described by 

where f{r,t) is the neutrino distributi on function at the position, r, and time, t, after the 



initial time, to, fiSwesty &: Myra 



20091 ). The diffusion coefficient, D, is related with the 
mean free path for the scattering, X, as D = cA/3. The parameter, /q, is the initial height 
of the packet. The value of to = d^/AD also gives the time scale of diffusion. The power 
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index, a, is related with the space dimension, A^^i, as a = No/'2- 

We first describe the spherical spreading of the Gaussian packet located at the center. 
The computation was done in 3D. This is a very crude approximation to what happens in 
the opaque region in the supernova core. We additionally examine non-radial diffusions 
of the Gaussian packet in the box geometry. Although this has no counterpart in the 
supern ova explosion, it is a common benchma rk test for neutrino (or radiation) transfer 



codes (ISwesty &: Myra 



2009 : 



Stone et al. 



I992I). 



For the first test we place the Gaussian packet at the center of the sphere with the 
radial coordinate extending from r = 0tor = 5x 10^ cm and the polar and azimuthal 
coordinates from ^ = to vr and from = to 27r. The initial width of the Gaussian packet 
is set to (io = 1 X 10^ cm. We consider only the isotropic scattering with a mean free path of 
A = 10^ cm. This corresponds to a diffusion time scale of 2.5 x 10~^ sec. We deploy A^^ =20, 
40, 80 and 160 radial grid points and Ng = 18 polar grid points and A^^ = 36 azimuthal grid 
points as spatial grids. For the momentum space, on the other hand, we employ = 2 for 
energy grid points and Nq^ = 12, N^^ = 12 for angle grid points. 

We show in Fig. [2] the radial profiles of neutrino density at the initial (t = s) and final 
{t = 1 X 10~^ s) time steps (10^ steps) for the case of A^^ = 80. The radial profile of number 
flux is also shown for the final time step. We compare the numerical results (cross symbols) 
with the analytic solutions (solid lines). It is evident that they agree well with each other. 
In Figure [31 we give relative deviations of the numerical results from the analytic solutions 
averaged over all the radial grid points for different grid sizes. As the number of radial grid 
points increases, the deviation is reduced as as predicted for the second order central 
differencing scheme. 

The second test for non-radial diffusion is done in 2D and 3D. In 2D, we place the 
Gaussian packet in the square that has a side of 10^ cm and is located at r = 10^ cm. In this 
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small area, coordinates are approximately Cartesian. We consider again only the isotropic 
scattering with the same mean free path of A = 10^ cm. We take A^^ = 100, Ng = 96 for the 
spatial grid and Ng^ = 12, A^^,^ = 12 for the angle grid points in the momentum space. We 
deploy A^e = 4 energy grid points despite the calculation is independent of neutrino energy 
in this test. We show in Fig. |4]the neutrino density at an early time {t = 5.0 x 10^^ s) 
and the final time (t = 5.0 x 10~^ s). The polar axis is denoted as Z (Z = rcos^) and 
the distance from the polar axis as R (i? = rsin^) in the plot. We compare the numerical 
result with the analytic solution along Z = —8.1 x 10^ cm (near the center of square) in 
Fig. 13 The deviation of the numerical result from the analytic solution is less than 0.5 % 
in the central region after 10*^ time-steps (At = 5 x 10^^ s). 

In order to see the convergence of the numerical results quantitatively, we repeat 
the same computation with different resolutions: (A^^, ^e) = (25,24), (50,48), (100,96) and 
(150,144); for the momentum space the grid points are essentially the same Ng^ = 12, 
N^^ = 12 except the number of energy grid points is reduced to A^e = 2 to save the memory 
size. In Figure El we show the relative deviations of the neutrino density averaged over the 
radial grid points along Z ~ 0. The deviation decreases quadratically as the number of grid 
points increases as expected for the central differencing scheme. 

The test is also done in 3D. We put a cubic box with a side of 10^ cm at a radial 
position of r = 10^ cm. The box is small enough to regard the patch of spherical coordinates 
as Cartesian. The mean free path of the isotropic scattering is again A = 10'^ cm. The 
numbers of grid points employed in this test are A^^ = 50, Ng = 48, A"^ = 48 for the spatial 
grid and Ng^ = 12, N^^ = 12, A^^ = 4 for the momentum space. We show in Fig. [7] the 
neutrino density on the plane x = 10^ cm for the initial (0 s) and final {t = 3.3 x 10~^ s) 
times as a contour map. The isosurfaces are also shown in color. We compare the numerical 
result with the analytic formula, Eq. fl22|) . The relative deviation is shown in Fig. |8l 
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The final distribution obtained by the numerical computation agrees very well with the 
analytic one. The deviation is typically less than 5 % in the central region for the case with 
= 50, Ng = 48, N^ = 48att = 5.0 x 10"^ s after 10^ time-steps (At = 5 x 10"* s). Large 
deviations near the edge of the box are mainly due to the fixed boundary condition we use 
for this test. 

In order to see the convergence of numerical solutions, we perform another test 
computation with a higher spatial resolution, A^,. = 75, Ng = 72, A^^ = 72, and a smaller 
time step of At = 2 x 10^*^ s. The results for the different grid sizes are compared at 
the same timing t = 2.0 x 10^^ s. The accuracy is improved by a factor of ~2 over 
the whole region. Unfortunately, further computations with even higher resolutions are 
impossible owing to the limitation of available memory space (See §A.5p and the power 
index of convergence cannot be determined at present. It is encouraging, however, that the 
numerical solution is indeed converging to the exact solution in a way that is consistent 
with the ID and 2D counterparts. We note incidentally that the total number of neutrinos 
in the box is conserved within the accumulation of round-off errors as expected for the 
conservative scheme. 

4.1.2. Free Streaming Limit 

We next examine the free streaming limit. We switch off all neutrino reactions in the 
following. As the simplest test, we compute a ID advection of neutrinos. We utilize a 
small radial interval located at a large distance. We set up a step-like distribution initially 
(top panel of Fig. [9]). All neutrinos have = 0.93247, which corresponds to the forward 
grid point for the polar coordinate of neutrino momentum whereas they are uniformly 
distributed in the azimuthal direction. The numbers of mesh points are A^^. = 100, Ng = 3, 
A'^ = 3 for the spatial grid and Ng^ = 6, N^^ = Q, N^ = A for the momentum space. This 
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means that we employ the full 3D code to solve the ID advection in space. The time step 
of 1.0 X 10~^ s is adopted to follow the evolution over 2.0 x 10^^ s. 

As shown in Fig. [9l the step in the neutrino distribution propagates radially at the 
projected light speed, ^^c- The edge of the step smears out gradually owing to numerical 
diffusions as it propagates. The upwind scheme is adopted for the current test (See \A.3\} . 

We next confirm that this numerical diffusion is reduced when we deploy finer spatial 
grids and take smaller time steps. We show in Fig. [10] the neutrino distributions at the 
final time (2.0 x 10~^ s, corresponding to the bottom panel in Fig. [9]) for different spatial 
and temporal resolutions. In the upper panel, we show the cases for different time steps, 
At=10^^ s, 10^^ s and 10^^ s (corresponding to the Courant numbers, 3, 0.3 and 0.03, 
respectively) for the same number of radial grid points of A^j.=100. The smearing becomes 
apparently smaller when we reduce the time step from 10"^ s to 10^^ s. The improvement 
is not so drastic when we take 10~^ s instead of 10^^ s. The further reduction of numerical 
diffusions is obtained by improving the spatial resolution. When we increase the number of 
radial grid points to A^r=1000 from Nr=100, the smearing of the step becomes even smaller 
as we can see in the lower panel. We note that the time step is simultaneously reduced by 
a factor of 10 in this computation so that the Courant number remains to be 0.3. 

The smearing of sharp edges in the dis tribution function as observed above is 



rather common to ra diation transfer codes (iStone et al. 



1992 



Turner fc Stone 



2001 



Swesty fc Myra 



20091 ). The performance of new code is not so bad as it seems, being 
comparable to those obtained by other popular transfer codes. For example, the reduction of 
numerical diffusion with the spatial and temporal r e soluti ons in this test is similar to what 



was shown in §4.3.3 with Fig. 16 of 



Swesty fc Myral (120091 ). No oscillation (or overshooting) 



around the edge in our results is found as expected for the implicit time-differencing 
we adopt in our code and is consistent with the results reported in §6.1 with Fig. 7 of 
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Stone et al.l (119921 ) 



Although the simple advection of the step-like distribution is done as a standard test 
here, we stress that it is too stringent from a view point of core-collapse simulations, since 
we do not obtain such a sharp edge in the neutrino distribution function as we will see 
in §5.11 for example. It is true that the smearing of the forward-peaked distribution in 
the optically thin region at large radii in the core is a concern. This is a problem more 
closely connected with the numbers of angular grid points in the momentum space as we 
will discuss in §4.2.11 rather than that of the advection scheme, though. We note also that 
it is the neutrino transfer and its influences to the hydrodynamics of material below the 
stalled shock wave (~200 km) that we would like to address with the new code. Then the 
forward peak of the neutrino angular distribution is not so appreciable and the numerical 
smearing will be less problematic. The neutrino luminosities and spectra at much larger 
radii are certainly important particularly from an observational point of view and will be 
addressed quantitatively by our new code on supercomputers of the next generation. The 
current advection scheme is admittedly diffusive and is adopted just as a first step in the 



3D neutrino transfer, which itself has beco me possible only recent 



Exploration of a more sophisticated scheme (jStone fc Mihalas 
future issue. 



y and is in its infancy. 



19921 . for example) is hence a 



As a final test, we present 3D computations of a searchlight beam (jStone et al 



19921 ). We note again that this test is also too severe for the code intended for supernova 
simulations in spite of its popularity as a benchmark test for radiative transfer codes. 
Indeed, there is no hot spot in the supernova core unlike in the sun which would require 
high angular resolutions. The main purpose of this test is to examine whether the code can 
give a correct propagation velocity of neutrinos, that is, the light speed in the genuinely 
3D setting. We inject the neutrino beam with fi^, = —0.90412 and (p^, = 6.1771 radian at 
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a certain point on the boundary. The numbers of mesh points employed for this test are 
Nr = 50, Ne = 48, = 48 for the spatial grid and Ne, = 12, N^^ = 12, A^'^ = 4 for the 
momentum space. 

In Figure [TTl we show the resultant neutrino densities in a 3D box with a side of 
10® cm at the time of ~lxlO^^ s when they become almost steady. We can see that the 
beam propagates along the designated direction, while the beam becomes broader as it 
propagates because of numerical diffusions. It is to be noted that the diffusion is mainly 
due to the small numbers of angular grid points in the momentum space rather than due to 
the relatively coarse spatial resolution, (see the employed grid in the 3D box presented in 
the upper panel). This is different from the situation for the ID advection discussed earlier. 
We also find that the beam shape is deformed as it propagates owing to the non-uniform 
angular grid in the momentum space along the beam. These results indicate that it is 
desirable in principle to deploy an angular mesh in the naomentuni spac e that is as fine as 



possible and covers the whole sohd angle uniformly (See 



Ott et al 



2008 



for example). 



4.2. Collision Term 

In order to examine the collision term for the neutrino reactions in supernova cores, we 
investigate test cases with realistic profiles using the neutrino reaction rates described in 
§3.21 We concentrate here on basic tests under spherical symmetry. We will explore further 
the neutrino transfer inside the supernova core in ^ 

As a typical situa tion, we utilize supernova profiles from the spherical calculation 



( ISumiyoshi et al. 



20051 ). We take the profiles of density, temperature and electron fraction 
at 100 ms after the core bounce in the collapse of the 15M0 star (See also §5.ip . We adopt 
the radial grid points {Nr = 200) for r = ~ 1.4 x 10^ km from the original profile. We 
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cover the first octant of tlie spliere witfi tlie angular grid, tliougli we treat a splierical profile. 



The neutrino energy and angular grid s are determinec 
spherical simulations of supernovae by ISumiyoshi et al 



jy fo llowing the setting in the 



t00& ). The energy grid (A^"^ = 14) 



is placed logarithmically from 0.9 MeV to 300 MeV with a fine grid for high energy tails. 



4-2.1. Stationary Case 

We examine stationary cases through comparisons with the formal solutions of neutrino 
transfer. For this purpose, we take into account absorptions and emissions in the collision 
term and switch off scatterings. We check the neutrino propagation with the neutrino 
reactions in dense matter along a certain ray. This covers the intermediate regime between 
the opaque and transparent regimes. Therefore, this is a complementary check to the tests 
in the diffusion and free-streaming limits. 

We evaluate the neu trino distribution at a certain location by the formal solution 



f Mihalas fc Mihalas 



19991 ). The neutrino distribution, /(s), at the position, s, along the ray 



follows the transfer equation given by, 

-^/(^) = -x(^)/(^) + ^(^), (23) 

where x(s) and r]{s) are the opacity and the emissivity, respectively. The path length, s, is 
measured backward (the opposite to the neutrino direction n) from the position (s = at r 
to obtain /(r)) to the outer boundary {s = Sbc at r = Rout)- The neutrino distribution at 
the position, r, is obtained by the integral, 

f{s)= r e-^^'\{s')ds\ (24) 
Jo 

along the ray designated by and (j)^. The optical depth, r(s), along the ray is given by, 

r{s)= !\{s')ds'- (25) 
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We assume in deriving Eq. the incoming neutrino flux is zero at the outer boundary. 

We integrate numerically the above equations using the rates of emission and absorption 
on the grid points along the ray of neutrino propagation. At the same time, we perform 
computations of the neutrino transfer by solving the Boltzmann equation in 3D for a 
sufficiently long period to obtain a steady state solution. We fix the spherical profile of 
the supernova core at the post-bounce mentioned above by the spatial grid with Nr = 200, 
Ne = 3 and A^^ = 3. We set the neutrino angle grid with Ng^ = 6, 12, 24 and N^^ = 6. We 
treat here the electron-type neutrino only. 

We compare the energy spectra obtained by the two methods in Fig. [12] (a). We plot 
the spectra at the radial position of r =98.4 km for the neutrino direction with fi,y = 0.99519 
(the most forward grid point) for Nq^ = 24. The distributions accord well with each other 
for the wide range of neutrino energy, though there is a slight deviation at low and high 
energies. In Fig. [12] (b), we show relative errors of the distributions by the computation 
with respect to the formal solutions. For the case of Ng^ = 24, errors are less than 30% for 
neutrino energies, < 100 MeV. Errors become large for high energy tails above 200 MeV 
since tails of distributions become fairly small below 10^^°. As the energy goes below ~3 
MeV, errors increase beyond 10%. This is because the matter becomes transparent for low 
energy neutrinos with small cross sections. A high angular resolution is necessary to follow 
the propagation of neutrinos in the transparent situation. For the medium range of energy, 
which is most important in the supernova study, the transition from the opaque central 
core to the transparent outer layers is well described by the computation. 

We next examine the angular resolution by changing Ng^. In Fig. [12] (b), we plot 
also the cases of Ng^ = 6 and 12 for comparison. The neutrino directions at the most 
forward grid points are fi^, = 0.93247 and 0.98156, respectively, in these cases. As Ng^ 
increases, errors become small in principle, showing the improvement of the angular 
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resolution. Although Ng^ = 24 (or larger) is preferable for the precise evaluation of 
forward-peaked distributions, Nq^ = 6 is sufficient to obtain errors less than 40%. Note that 
No.. = 6 is the minimally proper size for the supernova study as checked by the detailed 



tests ( Yamada et al. 



supernovae so far (ISumiyoshi et al. 



19991) and is adopte d for the spherical calculations of core-collapse 



2005 



for example). 



The angular resolution to describe the peaked distribution is the intrinsic problem 
of the Sn method in the field of radiation hydrodynamics. In the current problem of 
core-collapse supernovae, though, it is rather important to describe the phenomena around 
~200 km, where the shock wave is hovering, with the neutrino emission from the central 
core (~ 50 km). Therefore, the resolution for the angle factor of ~0.25 may be sufficient 
to describe the phenomena such as neutrino heating, at least for the first trial in 3D 
simulations. This is totally different from the situation in the solar physics, for example, 
where a small hot spot may be crucial for the radiation in outer layers at very large 
distance. Under the reasoning described above, we adopt Ng^ = 6 for the following study of 
supernova cores within the currently available computing resources. 



4-2.2. Time Evolution toward Equilibrium 

We examine the time evolution of the neutrino distribution toward the equilibrium 
determined by the condition of dense matter. We check the time scale to reach the 
equilibrium and the detailed balance under the chemical and thermal equilibrium through 
the neutrino reactions. We follow the evolution from an initial neutrino distribution with 
small values (10~^ times the equilibrium value) to the equilibrium state. In the static 
background, the time evolution of the neutrino distribution at a certain energy grid point 
by the absorption and emission can be analytically expressed as / = (/o — /fz))^"'^ + fpD, 
where /o and fpo are the initial and equilibrium (Fermi-Dirac) values, respectively. The 
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time scale, r, is given by r = ^ with the effective mean free path, A. The effective mean free 
path is defined here by the inverse of the true absorption coefficient, k,. The true absorption 
coefficient is given by the sum of contributions 

Kabs = Rabs{e,m^ + e-P^'-^^^] , (26) 

for electron captures and neutrino absorptions and 

^ j dn' R,air-emUe, s', U') [c'^^^^^' ^^0 + [1 " /(^', ^')]\ , (27) 



for the pair processes. Since we consider the iso-energetic scattering, the neutrino scattering 
contributes only to the realization of isotropy. The scattering coefficient, (Tscatt, is given by 

-2 



^scatt = TTTIyS / Rscati^',^') , (2^ 



which is larger than the absorption coefficient in this case. The effective mean free path 
due to the scattering process is short enough to realize the isotropic distribution within a 
short time period. We utilize the profile of the supernova core at 100 msec after the bounce 
as a background and choose a central grid point for the test. The neutrino reactions listed 
in ^3.21 are all included. 



In Figure [T31 we show the time evolution of the neutrino populations of z/g for the 
energy grid points at Ei,= 34.0 and 129 MeV at the center of the supernova core. The 
density, temperature and neutrino chemical potential are 3.15x10^^ g/cm^, 13.4 MeV and 
158 MeV, respectively. The time evolution in the computation is well in accord with the 
analytic solution. Time steps can be very long (over 1 sec) once the distribution reaches 
the chemical equilibrium owing to the implicit treatment. We show also the energy spectra 
in Fig. [131 The spectrum evolves toward the equilibrium and reaches the Fermi-Dirac 
distribution determined by the temperature and chemical potential. We note that the angle 
distribution becomes isotropic through scatterings during the evolution. 
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We check also that the approximation for the pair process, which is taken for the 
hnearization, is appropriate in the reahstic profile of supernova core. In the central part of 
core, where the temperature is high, the neutrino distribution approaches soon the thermal 
equilibrium as shown above. The blocking factor in the reaction rate for the pair process 
can be hence expressed by the thermal distribution or the distribution at the previous 
step. The effective mean free paths due to the pair processes, Eqs. (ITS]) and (IT^ . in the 
new code are compared with those from the spherical calculation. They accord very well 
with each other once the thermal equilibrium is maintained after a short period (See §5.ip . 
This time scale is much shorter than the hydrodynamical time scale (~ms), therefore, the 
approximation can be safely used in the dynamical simulations of core-collapse supernovae. 



5. Applications 

We investigate here the performance of our code for realistic profiles taken from 
supernova cores. We first employ two spherically symmetric core profiles, one during the 
collapse and the other for the post-bounce stage. In 2D and 3D cases, we deform by hand an 
originally spherically symmetric core profile rather arbitrarily and investigate the neutrino 
transfer in the non-spherical settings. We demonstrate that the new code can describe such 
features as fluxes and Eddington tensors in a qualitatively correct way both in 2D and 3D. 
We also use the formal solutions for more quantitative assessments. 



5.1. ID Configurations 



We examine first the neutrino transfer under spherical symmetry through comparisons 
with the numerical res ults by the neutrino-radiation hydrodynamics in general relativity 
( ISumiyoshi et al.l 120051 ). Adopting the profiles of supernova cores, we follow the time 
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evolution of the neutrino distributions from small initial values until they reach a steady 
state. We treat the three neutrino species, z/e, Ve and and implement the neutrino 
reactions listed in §3.21 



Figure [H] sho ws the adopted profi. " 



of the 15 Mq star ( Woosley fc Weaver 



es of the su pernova core in the gra- vitational collapse 



Sumiyoshi et al. 



(120051 ) ■ The snapshots 



19951) from 

are taken at the timings when the central density is 10^^ g/cm^ during the collapse and at 
100 ms after the core bounce. The former is an example of the situation during the collapse, 
where the neutrinos are trapped by the reactions with nuclei. The latter is a typical 
situation of the stalled shock wave after the bounce, where the neutrino heating takes place 
by the neutrino flux from the central core. We set the spatial grid with A^^ = 200, Nq = 5 
and = 5 in the first octant and take the radial grid points from the original profiles as 
done in §4.21 The neutrino angle and energy grids are set with Ng^ = 6, N^^ = 12 and 



14. 



We show the radial distributions of the number densities and fluxes for the three 
species of neutrinos in the left panel of Figs. [TS] and Uni for the profile during the collapse. 
We show the corresponding distributions for the profile after the bounce in the left panel 
of Figs. [T7| and [181 The degenerate neutrinos (ue) and thermal neutrinos (z/g and u^) at 
the central region are properly described with the tail of free-streaming fiuxes in the outer 
region. 

In Figures [T51 [ini [13 and [TSl we compare the current result s (cross symbols) with 



Sumiyoshi et al 



fl2005h . 



the numerical results by the spherical code (solid lines) based on 
Since we treat the steady state, we evaluate separately the neutrino distributions for the 
static background by solving the general relativistic Boltzmann equation under spherical 
symmetry (ID). Relative errors in the neutrino densities and fiuxes between the two 
evaluations are shown in the right panel of each figure. In general, the numerical results 
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accord very well with each other, while the errors of the density amount to a few tens of 
percent near the boundary. The errors of the fluxes are significant in the region where the 
errors of the density are appreciable. We found that the general relativistic treatment in 
the ID code does not affect the comparison. This is because the neutrino distributions are 
determined locally by the thermodynamical condition at the central region and the neutrino 
fluxes at the outer layers are not affected by the general relativity due to large radii. 

We show in Fig. [19] the mean free paths of neutrino reactions in the profile after the 
bounce described above. The mean free path we discuss hereafter is the effective mean free 
path defined by the inverse of the absorption and scattering coefficients in Eqs. (12^ . (I27|) 
and (l28l) for each process. We check the mean free paths by comparing with the numerical 
results obtained by the ID code. As shown in Fig. [19], the mean free paths by the 3D code 
accord with those in the ID evaluation. The relative errors between the two evaluations are 
within 2 x 10^^ except for the cases of the pair process and nucleon-nucleon bremsstrahlung 
to be discussed below. The mean free paths for the collapse phase have been also checked 
in the same way (not shown here in figure). 

We note that the mean free paths for the pair processes within the approximate 
expression accord very well with the full expression in the ID code in the central region 
where these reactions are important. We show in Fig. [20] the ratios of the mean free path 
by the 3D calculation to that by the ID code for the pair process and nucleon-nucleon 
bremsstrahlung. In the central region, the ratio is very close to 1 and its deviation is 
within 4 x 10^'^. Deviations in the outer layers appear due to the approximation of the 
reaction rate and the angle-averaged distribution of couterpart-neutrinos. Deviations for 
the nucleon-nucleon bremsstrahlung are rather small since we adopt the same isotropic rate 
in the both codes. Large deviations in the pair-process are due to the approximation of 
taking isotropic reaction rates. We note, however, that the deviations in the outer layers 
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are not crucial for the whole description since the material is in the transparent regime and 
contributions of these reactions are minor there. 

In Figure [21], we show the emissivities for the three neutrino species as a function of 
radius for the profile after the bounce. The emissivities are defined by 

f de f 

{e,Q), (29) 

and 

[ de £^ [ f de'e'"^ [ , ^ ^ _ ^ ^ 

(^emis—pair — J ^27i-)3 J J (27r)^ J ^ ^pair—emisi^ ■> ^] ^ 5^ ) [l /(^ )^ )]> ("^0) 

for neutrino-emissions (ecp, aecp and eca) and pair-processes (pap and nbr), respectively. 
We compare the emissivities with those obtained by the ID code. The emissivities in the 
two evaluations agree very well with each other within relative errors of 2 x 10"'^ in general. 
The largest errors arise in the transitional region between 20 km and 100 km but are within 
~ 1 X 10^^. The agreement is good even for the pair processes despite the usage of the 
isotropic rates. This is because the isotropic term of the reaction rates is dominant in 
Eq. (!30|) with isotropic distributions in the inner region or small distributions in the outer 
region. This situation is different from the case of the effective mean free path for the pair 
processes, where the angular distribution is important due to the exponential factor in Eq. 

m 



We examine the energy and angle moments of the neutrino distributions obtained in 
the 3D code through the comparison with the ID evaluation. We show in Fig. [22] the 
Eddington factors, k^^, and the flux factors, {fiu), as functions of radius in the post-bounce 
profile. The definition of various moments of the neutrino distributions is summarized in 
§A.4[ We obtain the correct limits of quantities in the opaque and transparent regimes. 
The Eddington factor and the flux factor are 1/3 and zero, respectively, in the central core, 
where the distributions are isotropic, and they approach ~1 (forward peaked) toward the 
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outer layers. The moments by the 3D code accord generally well with those by the ID 
evaluation. In Figure [221 we show differences of the moments by the 3D code from the ID 
evaluation. Deviations within ~0.05 appear around the transitional region and near the 
boundary. 

We show, in Fig. [2U the luminosities and average energies of neutrinos for the 
post-bounce profile with the corresponding quantities from the ID code for comparison. We 
plot here the average energy, (e), and the luminosity, L^, = Anr'^Fiy, defined in §A.4[ The 
behavior of the luminosities and average energies accords generally well with the results by 
the ID code. 

The numerical checks so far prove the new code with the microphysics of neutrinos in 
dense matter works properly with reliable accuracy in the spherical configurations. From 
the examination of the effective mean free paths and the emissivities, we judge that the 
approximate expressions adopted in the pair processes are efficient and sufficient for the 
numerical studies of supernova cores by the 3D code. 



5.2. 2D Configurations 



In order to demonstrate the ability of the new code in mult i- dimensional realistic 
settings, we first study 2D case, utilizing artificially deformed p rofiles of supernova core 



based on the ID core-collapse simulation f lSumiyoshi et al 



20051 ). For the given background 



profile, we obtain a steady state solution of neutrino transfer by following the time evolution 
for a sufficiently long time. We compare the result with the formal solutions as discussed in 



We utilize the same spherical profiles for the post-bounce stage as those studied in §5.1 
We modify the profiles of density, p(r), temperature, T(r), and electron fraction, Ye(r), 
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along each radial direction, depending on its polar angle, 6, by scaling the radius, r, as 

r = r(l — esin^), (31) 

where e is a parameter to specify the degree of deformation. As a result, we obtain oblate 
profiles for positive e's, which are crude approximations to rapidly rotating supernova cores. 
We show in Figs. and 121] some features of the profiles thus constructed for e = 0.4. 

We first examine the neutrino transfer through comparisons with the formal solutions. 
For this purpose, we treat the electron-type neutrino alone, taking into account only the 
emissions and absorptions. We set the spatial grid with A^^ = 200, Nq = 9 and = 9 in the 
first octant. We try four sets of angular meshes in the momentum space: {Ng^, A'<^^) = (6,6), 
(12,6), (24,6) and (12,12). The number of energy grid points is fixed to = 14. The 
relatively small grid sizes are mainly due to the limitation of available computing resources. 

The formal solution is obtained by integrating the Boltzmann equation along the 
representative paths shown in Fig. [221 These paths go through one of four grid points that 
are rather arbitrarily chosen and have a radius of r=98.4 or 198 km with an angle cosine 
of /i = 8.2 X 10^^ or /i = 0.98. The solid lines in the figure represent the paths with a 
polar angle cosine in the momentum space of fi^, = 0.99519, which corresponds to the most 
forward grid point for Ng^ = 24. For comparison, the long-dashed and dashed lines show 
the paths with fi^, = 0.93247 and fi^ = 0.98156, which correspond to the most forward grid 
points for Ng^^ = 6 and 12, respectively. 

Figures |27] and [28] show the comparison of the numerical results by the new code with 
the formal solutions. In Figure [271 "we present the energy spectra for the points close to 
the equator (square symbols in Fig. [261) . The numerical results agree with the formal 
solutions within relative errors of 20 % at r=98.4 km. In this case, the density at the 
grid point is rather high (~10^^ g/cm^) and the matter is already opaque to neutrinos 
there. Moreover, the path runs through the proto-neutron star. As a result, the neutrino 
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spectrum at the grid point is mainly formed by the contributions from relatively near-zones 
and the agreement of the numerical result with the formal solution is excellent even with 
the rather coarse grid adopted in this computation. At r=198 km, on the other hand, the 
agreement is not so good as shown in the right panel of Fig. [271 This is because the density 
at the grid point is much lower and the matter in the vicinity is transparent to neutrinos 
as well as because the path barely touches the periphery of the proto-neutron star. Then 
the neutrino distribution at the grid point is a superposition of the contributions from very 
far regions. In general, the larger the distance to the source is, the more anisotropic the 
neutrino distribution becomes and the more difficult it is to numerically reproduce it. It is 
important to see, however, that the relative error is still within ~50 % even at very low 
energies, where the opacities are the lowest. 

These features are common to the points near the pole (triangle symbols in Fig. f2E\i 
as seen in Fig. |28l the energy spectra obtained numerically agree with the formal solutions 
within ~50 % except for the high energy tail beyond ~100 MeV, where the neutrino 
populations are very small, < 10^^'^. The relative errors near the pole are somewhat 
larger in general than those near the equator. This is because of the oblate shape of the 
artificially deformed core in this computation. The matter densities at the two points near 
the pole are lower than those at the corresponding points near the equator. The matter is 
hence more transparent in the vicinity of the former, hence, the evaluation of the forward 
peak is further difficult. The paths near the pole, on the other hand, run into deeper inside 
the proto-neutron star. Although these effects are counteracting each other, the former 
turns out to be more important. 

The relative errors become smaller for Ne^ = 24 than those for Ne^ = 6, 12. The 
increase of N^p^ to 12, on the other hand, does not improve the accuracy very much for 
Ng^ = 12. This feature is common to all the paths considered here. Although we need 
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further systematic studies with finer grids, the above facts suggest that it is more important 
to employ a sufficiently large Nq^ in the momentum space. 

We proceed to the numerical results with the full set of neutrino reactions for the three 
neutrino species. We set the number of grid points to Nj. = 200, Nq = 9, = 9, Ne^ = 6 
and N^^ = 12. We follow the evolution for a sufficiently long time period (~10 ms) to 
obtain a steady state. We show in Fig. [29] the contour plots of the neutrino density on 
the meridian slice with a constant 0=0.44 radian. The neutrino distributions apparently 
reflect the deformed proflles of density and temperature. The electron-type neutrinos are 
mostly populated in the central region, while the electron-type anti-neutrinos and mu-type 
neutrinos mainly exist off-center, where the temperature is high. 

Figure [30] presents the radial proflles of the neutrino density and flux along the 
directions with fi = 8.2 x 10~^ (near the equator) and /i = 0.98 (near the pole). The 
peaks of the densities for electron-type neutrinos are located at the center with different 
widths, reflecting the deformation of the proto-neutron star. The peaks for electron-type 
anti-neutrinos and mu-type neutrinos, on the other hand, are located off-center at different 
radial positions corresponding to the temperature peaks. These neutrinos are mainly 
produced by the pair-processes in non-degenerate and positron-abundant environment. 
This is the reason why they are abundant off-center. We flnd that their populations agree 
very well with the local equilibrium distributions. The radial neutrino fluxes also reflect the 
deformed matter distributions. In fact, the radial fluxes near the pole are larger than those 
near the equator at the same radius, since the density gradient is steeper near the pole (See 

Fig. m- 

We show not only the radial neutrino flux but also the polar flux of electron-type 
anti- neutrinos in Fig. [HI] to elucidate multi-dimensional transfer. The radial flux is enhanced 
near the pole as just mentioned. The polar flux, on the other hand, is signiflcant in the 



-40 - 



middle region between the pole and the equatorial plane. Because of the deformation, the 
density gradient is directed to the pole in general and is greatest in the middle region, since 
the axial symmetry and equatorial symmetry imposed in this computation force the density 
gradients radially directed near the pole and equator. We stress that this is a feature that 
can be captured properly only by the multi-dimensional transfer computations such as done 
in the current study and not by the ray-by-ray approximation. We note that the polar flux 



is significant even 



jeyond 100 km, w here an approximation by the polar gradient of the 



neutrino pressure (iMiiller et al. 



20101) may break down. 



We examine the energy and angle moments of neutrino distributions in Figs. [5^ and 
The flux factors, (/ii^), are shown in the upper panel of Fig. [32] as functions of radius along 
the two directions discussed above. The flux factors change from at center to 1 at large 
radii for both directions. The transitional zone corresponds to the region with the optical 
depth ~1 and has a different radial location for each direction. The average energies, (e), 
are shown in the lower panel of Fig. |32] as functions of radius. The radial dependence of 
energies for the three species of neutrinos reflects the thermodynamical states (degenerate 
or not) in the deformed core. The profile of the average energy for electron- type neutrinos 
follows roughly the density profile. The average energy declines rapidly as the density 
(and hence the degeneracy parameter given by chemical potential divided by temperature) 
decreases with radius. The profiles for other types of neutrinos (electron-type anti-neutrinos 
and mu-type neutrinos), on the other hand, reflect the temperature profile, having a peak 
around 10-20 km. The radial profiles near the equator are shifted outwards from those near 
the pole just as the other quantities seen above. 

As another measure of the transition from the opaque to transparent regimes, we 
show the Eddington tensor as a function of radius along the same two directions in Fig. 
The diagonal elements of the Eddington tensor (rr-, 66- and 00-components) are 
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shown for mu-type neutrinos with three different energies. The diagonal elements are 
1/3 in the central region, where the matter is opaque and the neutrino distributions are 
isotropic. The rr-component increases with radius as the neutrino distributions become 
more forward-peaked, whereas 99- and 00-components decrease. All off-diagonal elements 
are found to be nearly zero at the central region and have small values also at large radii. 
The transition of the diagonal elements from 1/3 to 1 or occurs around the neutrino 
sphere, which has larger radii for higher neutrino energies. 

We remark that the energy dependence is not so strong for electron-type neutrinos and 
anti-neutrinos and the radial profiles of the Eddington tensor are more close to each other 
for the three neutrino energies (not shown in the figure). This is because the isotropy of the 
neutrino distributions is maintained up to ~100 km through the charged current process. 



Ottetal 



(120081). Although the density 



This behavior is consistent with the analysis by 
profile is rather more extended in our model compared with theirs, they indeed found the 
delayed decoupling near the equator in their rotating model as seen in the right panel of 
their Fig. 6. 

The radial distributions of the Eddington tensor also depends on the polar angle (see 
top and bottom panels of Fig. |33]). Because of smaller density scale heights near the pole 
in the deformed core, the transition from the opaque to transparent regimes occurs at 
deeper and narrower locations near the pole as seen in Fig. |25l These are just as expected 
intuitively for the oblate core and demonstrate that the new code works appropriately at 
least qualitatively for the 2D configurations considered here. 
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5.3. 3D Configurations 

We demonstrate here the performance of the new code in 3D realistic settings. We 
study the 3D neutrino transfer utihzing deformed profiles in a similar manner to those 
studied in §5.21 We modify the deformed profiles further by adding the dependence on 
azimuthal angle, 0, in the scaling as 

f = r(l — esin0sin^). (32) 

The resulting 3D profiles are deformed maximally with polar dependence in the yz-plane 
(0=7r/2), whereas they have no polar dependence in the zx-plane (0=0). These profiles are 
simple examples inspired by the as ymmetric shape of supernova cores in the 3 D standing 



accretion shock instability (SASI) (IBlondin &: Shaw 



2007 



Iwakami et al 



20081). 



We first check the neutrino transfer through comparisons with the formal solutions 
in the same way as in 2D. We treat the electron-type neutrino with the emissions and 
absorptions. We set the spatial grid with Nr = 200, Ng = 9 and A^,^ = 9 and the angle 
grid with Ng^ = 6 and N^^ = 6. Figure [31] shows the comparison between the steady and 
formal solutions. Four panels show the energy spectra at four locations with r=98.4 km 
in different directions. The left panels correspond to the spectra for two polar directions, 
fi = 8.2 X 10^^ (top) and 0.98 (bottom) on the meridian slice with 0=0.262 radian (near 
the zx-plane). The right panels correspond to the spectra for fi = 8.2 x 10^^ (top) and 0.98 
(bottom) on the slice with 0=1.309 radian (near the yz-plane). The agreement between the 
two solutions is generally good for the wide range of energy, having relative errors within 
~20 % at energies around 20 MeV. Larger errors at low and high energies arise due to the 
same reason as described in §4.2.11 The energy spectra in the right panels differ each other 
due to the deformation with polar dependence. The energy spectrum near the equator 
(right top) is harder than that near the pole (right bottom) due to different densities and 
temperatures at the two locations. The two energy spectra in the left panels are similar 
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each other due to a nearly spherical geometry of the background. 

Next, we examine the numerical results with the full set of neutrino reactions. We set 
the angle grid with Nq^ = 6 and A^^,^ = 12. We display in Figs. 135) and 136) the density 
and fluxes of electron-type anti-neutrinos on the slices with 0=0.436 and 1.309 radian, 
respectively, by color contour maps. Figure demonstrates the deformations of density 
and flux profiles near the yz-plane, having the background with strong polar- dependence. 
The deformed profile of neutrino density follows exactly the shape of the deformed profile of 
temperature. The radial flux is enhanced around the pole and the polar flux is appreciable 
in the middle region as discussed in §5.2[ The azimuthal flux is not significant due to a 
small gradient in the azimuthal direction near the yz-plane. In contrast, the density and 
radial flux near the zx-plane in Fig. [35] have less deformed shapes than those near the 
yz-plane. While the polar flux has a certain contribution in the middle, the azimuthal flux 
has a significant magnitude due to the azimuthal dependence of the background near the 
zx-plane. 

We stress that our code properly describes the polar and azimuthal transfer, reflecting 
the deformation of the 3D supernova core. This is the advantage of the 3D Boltzmann 
solver in 3D profiles. The non-radial transfer cannot be described correctly in the ray-by-ray 
approach. We found that the polar and azimuthal fluxes are appreciable as compared with 
the radial flux in a wide region as seen in Figs. [35] and [36] The non-radial fluxes are 
significant even at r >100 km and spread beyond the diffusion regime. In this case, the 
reliability of the flux limited diffusion approximation is doubtful. Our code is, therefore, 
a unique tool to study the neutrino transfer in 3D configurations and to examine the 
approximations used in the state-of-the-art supernova simulations. 

We examine the energy and angle moments of the neutrino distributions in Figure [23 
Contour plots show the angle moments, (/x^), and the average energies, (e), for electron-type 
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anti-neutrinos on the slices with = 0.436 (near the zx-plane: left) and (j) = 1.309 (near 
the yz-plane: right) radian. Deformation is strong near the yz-plane and weak near the 
zx-plane due to the polar and azimuthal dependence. The angle moments have an isotropic 
value of 1/3 at the central region as shown by the deformed shape of contour lines. They 
approach ~0.8 at ~500 km independent of the polar direction as the angle distributions 
approach a forward-peaked shape. The contour lines for the average energy have also a 
deformed shape, reflecting the deformed supernova core. The average energies at large radii 
(~1000 km, not shown in the figure) do not show any significant asymmetry. 

Figure |3S] shows the elements of the Eddington tensor as functions of radius along 
the two polar directions, /x = 8.2 x 10^^ (top) and /i = 0.98 (bottom), on the slice with 
= 1.309 radian (near the yz-plane). Six elements of the tensor are shown for electron-type 
anti-neutrinos at energy grid point of 34.0 MeV. The transitions from the isotropic values 
(1/3 for diagonal and for non-diagonal elements) occur further locations in the direction 
near the equator than those near the pole. This dependence on the polar direction is 
small near the zx-plane and the transitions occur similar locations (not shown here). The 
magnitude of non-diagonal elements is appreciable at outer layers having the azimuthal 
dependence in the 3D supernova core. The transition from the isotropic regime to the 
free-streaming regime is described well by the new code as in the ID and 2D cases. Our 
code is also useful to examine the Eddington tensors for general 3D profiles and to provide 
valuable information for the approaches using variable Eddington factors. 

6. Summary and Discussion 

We have developed the numerical code to calculate the neutrino transfer with 
multi-energy and multi-angle in three dimensions for the study of core-collapse supernovae. 
The numerical code solves the Boltzmann equations for the neutrino distributions by 
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the discrete-ordinate (S„) method. The time step is advanced by the fully imphcit 
differencing. The neutrino distribution is a function of time, three spatial coordinates, 
neutrino energy and two angles. We solve, therefore, the time evolution of the distributions 
in six phase-space dimensions. An essential set of neutrino reactions (emission, absorption, 
scattering and pair processes) is implemented in the collision term of the Boltzmann 
equation. The Boltzmann equations are formulated in the inertial frame as the basis of 
our developments. In the current study for the static background, we drop the velocity 
dependent terms in the collision term. The collision term of the pair-processes is linearized 
by assuming the neutrino distributions for the counterpart of the neutrino pair. The set 
of equations for the discretized form of the Boltzmann equations is formulated as a linear 
equation for the vector of neutrino distributions. The large sparse matrix with block 
diagonal matrices is solved by an iterative method. 

We have performed the numerical simulations by the 3D Boltzmann solver to vahdate 
the numerical code and to apply it to the realistic profiles of supernova core. The numerical 
code of the 3D Boltzmann solver describes correctly the neutrino transfer in the diffusion 
and free streaming limits. The collision term by the neutrino reactions has been checked 
by the formal solutions as well as the analytic evolution toward the equilibrium. The 
pair processes are well described in an approximate expression for the linearization of the 
collision term. We have examined the neutrino transfer in the profiles of the supernova 
cores before and after the bounce using the numerical results from the previous simulations 
under spherical symmetry. The 3D calculations using the spherical configuration show the 
good agreement in neutrino distributions with the spherical calculations. 

We have demonstrated that the ability of the 3D numerical code to solve the neutrino 
transfer for artificially deformed supernova cores in 2D and 3D. The new code describes 
properly the multi-dimensional feature such as lateral and azimuthal fluxes, which cannot 
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be handled in approximate methods. The moments of energy and angle of neutrino 
distributions in 3D can be revealed by the new code and will be helpful to gauge simplified 
transfer approaches. 

Toward the full calculation of core-collapse supernovae by the neutrino-radiation 
hydrodynamics, there are necessary developments on top of the current 3D code. First 
of all, we need to link the 3D Boltzmann solver with the numerical code for the 3D 
hydrodynamics. This development for the neutrino-radiation hydrodynamics is now under 
way. In order to describe the neutrino transfer in hydrodynamics, it is necessary to take 
into account the relativistic effects {y/c terms) such as the Doppler shifts of neutrino energy 
and the aberration of neutrino angles. We plan to evaluate the collision term by considering 
the Lorentz transformation of the neutrino distributions between the inertial and comoving 
frames. This causes changes of neutrino energies by scatterings in the inertial frame even for 
the iso-energetic scattering in the comoving frame. Therefore, the consideration of energy 
change in the collision term becomes mandatory. Studies by the developments mentioned 
above will be reported as the following articles after the current study in series. 

Further generalization of the 3D numerical code would require supercomputers 
in Exa-scale. In principle, one has to include the energy- changing reactions such as 
electron-neutrino scatterings in the collision term. Considering the energy coupling makes 
the size of block matrix larger and the computation formidable. The non-linear treatment 
of the pair processes requires iterations of the whole solution. The general relativistic 
treatment of the whole framework makes the computation further challenging. Therefore, 
the full 3D simulations with the extensions mentioned above awaits the next generation 
supercomputers to perform. The current study proves that the calculation of the neutrino 
transfer in 3D is now feasible, marking the first step toward this direction. 
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A. Appendices 
A.l. Neutrino Direction 

We define the unit vectors at the position (r, 6, 0) as 



e^ = sin 9 cos 01 + sin 9 sin j + cos 9 k, 
eg = cos 9 cos i + cos 9 sin j — sin 9 k. 



—sin i + cos j. 



(Al) 
(A2) 
(A3) 



19731 ). The neutrino direction in 



The reversed transformation can be found in (iPomraningI 
the inertial frame is defined by two angles, 9^ and 0,^, with respect to these unit vectors. 
The polar angle between the neutrino direction, n™, and the radial coordinate, e,,, is 
denoted by 9^. The azimuthal angle around the radial coordinate from the e^-direction is 
denoted by (j)y. The three components of n*" are given by 



nT = cos 9„ 



n^Q = sin 9 1, cos (p^, 



nT = sin 9u sin 0,^. 



(A4) 

(A5) 
(A6) 



They can be expressed in terms of the angle variable as 



(A7) 
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[l-lil)^ COS 



sm 



(A8) 
(A9) 



A. 2. Computational Grids 

We arrange the radial, polar and azimuthal coordinates by N^., Nq and grid points. 
We cover the neutrino energy and angles by A^^, Ng^ and N^^ grid points. Hereafter, we use 
lowercase subscripts (e.g. i, j) for the mesh centers and uppercase subscripts (e.g. /, J) to 
index the location. We assign the integer index from 1 to A^ for the cell centers and the 
integer index from to A^^ for the cell interfaces, where A^^ is the number of grid points for 



Jth interfaces (e.g. 



Ng anc 




Yamada 


1997) 



The radial grid points at the cell interface are denoted by r/, which covers from r/=o to 
Ti=Nr- We define the radial grid points at the cell center by 



(AlO) 



We determine the grid points for the neutrino angle, Hu = cos 9,^, by the points and weights 
obtained from the formula of the Gaussian quadrature. The fi^-ghd points for the cell 
center, fii^j, are given by the quadrature points. The fiu-grid points for the cell interface are 
determined by 

l^iyj = l^-yj-i + d/J-uj , (All) 



where dfi^, are the quadrature weights for the cell containing /Xj, •. Starting with /i, 



we can determine the /ii^-grid points up to fiuj=Ng = 1. The angle factor at the cell 
interface, which is necessary for the evaluation of the /i,y-advection, is obtained by 



'1 - fiu'^)j-i - 2fi^Jfi^ 



(A12) 
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The angle factor at the boundary is set as (1 — iJiJ^)j=o — and those at other grid points 
are determined accordingly up to (1 — iJii?)j=Ne^ = 0- We use indices i and / for spatial 
coordinates and j and J for neutrino angle coordinates, hereafter. 

We arrange the grid points for the polar coordinate, 9, in terms of = cos 9 in the 
same manner as the /Xi,-grid points described above. We set the grid for the cell center, /Xj, 
and the corresponding weights, d/ii, from the Gaussian quadrature points. The grid points 
for the cell interface are determined by 

III = + d/ii , (A13) 

to cover from /x/^o to iJ>i=Ng- The angle factor at the cell interface for the polar advection is 

given by 

(1 - = (1 - - T-^^dn, , (A14) 

starting from the value at the boundary as (1 — fi'^)^ j^q = (1 — fii=o'^)^. Similarly, the 
01/- grid points at the cell center are set by the Gaussian quadrature with and dcp^i- The 
grid points at the cell interface are determined by 

01^7 = <l>i^j-i + d<l>uj , (A15) 

from (f)uj=o — to — n. The angle factor at the cell interface used for the advection 
is given by 

(sin0^)j = (sin0j,)j_i + cos^j,^- d(j)^j . (A16) 

starting with (sin0,^)j=o = 0. For three dimensional calculations, we need to cover from 
to 27r. We repeat the above procedure from = tt to = 27r after dividing the whole 
angle range into two parts. This is to avoid the 0,^-grid being asymmetric around the radial 
direction. 

The 0-grid is given in a similar manner to 0,/-grid. The 0-grid points for the cell center 
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are denoted by 0i and those for the cell interface are given by 

0/ = (Pi.i + . (A17) 
The points, (pi, and weights, dcp^, are taken from the Gaussian quadrature points. 
Finally, the energy grid points at the cell center are defined by 

Sk = {skSk-iY^^ , (A18) 

from the energy grid points at the cell boundary, Ek, covering the energy range from eK=Q 
to eK=Ne- The energy integration is done by the cell volume for the energy phase space, 
e^de = (i(£:^/3), as 

^(4/3) = (4-4-i)/3. (A19) 



A. 3. Numerical Scheme for Advection 



We describe a discretized form for the advection terms in our s cheme . The current 
approach is an extension of the method by iMezzacappa &: BruennI (jl993[ ) adopting an 



interpolation between first- and second-order finite difference representation of the neutrino 
advection. We pay attention to the evaluation of quantities at the cell interfaces for a wide 
range of opacities. We assume the evaluation is done at the time step n -t- 1, although we 
drop here the index. 

First, we explain the advection terms for r and ^i, which are the two basic variables 
used under spherical symmetry. The radial advection is expressed by 



d 



9(rV3) 



ir] fi-rUfi-i), 



(A20) 



i-i 



where and // are the neutrino distributions at the cell interfaces for /j. The value of 
f^ujfi at the cell boundary is evaluated by 

/X.,// = ^^^^^^^{(1 - + /3,/,+i} + + (1 - (A21) 



- 52 - 



for both inward and outward directions of neutrinos, depending on the sign of fi^j- For 
example, the expression for outward neutrinos (/x^^ > 0) becomes 

f^u.fi = + (1 - , (A22) 

where /3/ is a weighting factor to bridge the upwind and central differencing. We note 

that the expression, // = /«, is obtained for the upwind differencing by setting /?/ = 1, 

whereas the expression, // = (/j + /j+i)/2, is obtained for the central differencing by setting 

/3/ = 1/2. The value of /3/ is determined by a smooth function to connect the diffusion 

= 1/2) and free-streaming = 1) regimes. We employ the following expression 

1 aArj/Xj 



Pi 



(A23) 



based on the formula by 



?1 +aArj/Aj 

Mezzacappa fc Bruenru (119931 ). Here the interval of radial grid 



points is defined by Arj = rj+i — and the average of mean free paths is defined by 

= (Aj^\ + \^^)/2. An adjustable factor, a, is set to be 100 to ensure (3j = 1/2 for the 
diffusion regime even when Arj ~ A/. 

The advection term for fi^, = cos 6^ is given by 



1 d 



r dn. 



3 r] — r 



i-i 



[(1 - f^u')jfj - (1 - fiu')j-ifj-i] , (A24) 



2 r| - r]_-^ dn^j 

where we take the upwind differencing and set simply fj = fj. We use the angle factor 
given by Eq. flA12p . The relation of Eq. f lA12p must be used to guarantee the steady state 
for infinite homogeneous medium (constant /«) in the neutrino transfer under spherical 
symmetry by adopting the ad vection terms with Eqs. (1A20I) and (1A24I) as described in 



Mezzacappa fc Bruenn 



fll993h . 



Next, we describe the advection terms for 6 and 0^, which are essential for the 
description under axial symmetry. We rewrite the 6'-advection term by using the variable n 
and discretize the term by 

^/l- cos (l)ud_^^.^ _ [ fJ-l 

'de 



rsin 6 



[sin Of) 



cos 



-I 
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2 T'^j 



^ ''l-/^')l//.-(l-/^')Li//«-i 



«9 



(A25) 



We use the angle factor given by Eq. (IA14I) . The factor, (1 — ix^'^J^cos(f)uj^, determines the 
direction of advection and the evaluation of fig at the cell interface. Depending on the sign 
of cos0,^, the value at the interface is given by 



cos (f)^ J Jig 



COS0,,^ + |COS0,,.J 



COS( 



COS( 



(A26) 

This expression is adopted to connect the diffusion and free-streaming limits smoothly by 
the factor, Pig, in the same way as the one used in the radial advection. The form of Pig is 
given in a similar form as Eq. flA23|) by replacing the width of grids and the average mean 
free path to the ones in polar direction. 



The advection term for (p^y is given by 
yT^^T^cos 6 d 



sin 6 d(j)i. 



2 ij^^ 

Ir 



2 \i 



(sin 0^/) 



d 



(sin (p^f) 



1 



[(sin - (sin , (A27) 



using the variable, fi. The angle factor is given by Eq. ( IA16|) . The sign of /Xjg(sin0^)j^ 
determines the direction of advection. We take the upwind differencing for the evaluation 
of fj^ at the cell interface through the relation as 



2 

/ii,(sin</)^)j^ - |/i,,(sin0^)j^ 



(A28) 



In the same manner as the spherical case, the combination of Eqs. (lA14p and (1A16P 
guarantees the steady state for infinite homogeneous medium in the neutrino transfer under 
axial symmetry when the ad v ection terms of 9 and (py are expressed by Eqs. flA25P and 



flA27p as described in 



Castoil mm 
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Finally, we describe the advection term for 0, which is necessary for three dimensional 
calculations. The advection term is expressed in terms of n as 

a/I ~ 1^1 9f 



3 rf 



Ir~l 



2 r? 



Ir-1 



rsin 9 



(1-/^4) 



1 



sin 



df 



sin 



r-\/l — /x^ 
1 



The evaluation of fi^ is made by 



sm I 



sin 



sm( 



+ 



2 

|sin( 



{l^Ufis + (1 ~ /^-f,A)/u+l} 



(A29) 



(A30) 



depending on the sign of sin0,^j^. The form of f3i^ is given by a smooth function, which 
is similar to /3i^ and ^j^. This form of 0-advection fulfills the steady state in infinite 
homogenous matter and the advection vanishes when /j is constant. 



A. 4. Moments of Energy and Angle 

We define the moments of energy and angle of the neutrino distributions. All 
evaluations here are done using the quantities in the inertial frame, though we dropped 
superscripts for the compactness. First of all, the neutrino density is evaluated by 

f 7^ ! dnf{e,Q). (A31) 



J (27r)3 

The first moment of energy is defined by 



{e) = f^, (A32) 

77/1/ 



where the energy density of neutrinos is given by 

E.^ j ^ j dnef{e,n). (A33) 



-55 - 



The first moment of angle is defined by 



= ^, (A34) 



where the radial number flux is given by 



/:= / / dnnrf{s,n). (A35) 



(27r)3 

The polar and azimuthal fluxes are obtained by 

= dnnef{e,n), (A36) 

The radial luminosity of neutrinos is defined by 

U = '^nr^F^, (A38) 

where the energy flux is given by 

F. = J ^ J dn snrf{s,n). (A39) 

The second moment of angle is defined by 

= - dn nrnrf{e,n). (A40) 



n,J (27r)3 

We evaluate the elements of the Eddington tensor as defined by 



where the elements of the pressure tensor are defined by 



pij 

k'' = ^, (A41) 



Pi' ^ j^.jdn en,njfie,n). (A42) 



The subscripts, i and j, denote one of the three components of r, 9 and (p. The diagonal 
elements in the spherical coordinate are A;^^, k^^ and k'^'^. The non-diagonal elements, 
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f^re^ f^rcj) f^eip^ j-j^g^y non-zero in 3D calculations, while they are zero under spherical 
symmetry. For the treatment with multi-energy group, we utilize the integrand of Eqs. 
( IA33P and flA42p with the angle average. We evaluate the Eddington tensor for each energy 
zone by 

dropping off the common factor of the energy phase space. 



A. 5. Computing Size 

We describe briefly the size of memory and computational load necessary for the 
current simulations as well as larger simulations in future. The typical size of memory 
requirement for the numerical calculations in §5.21 and §5.31 is ~30GB in the case of 
Nr = 200, Ne = 9 and = 9 with A^'^ = 14, Ng^ = 6 and N^^ = 12. This includes the 
matrices for the equations and the vectors for the neutrino distributions for three species. 
We need ~130 MB to store the neutrino distribution for each species. It takes ~100 sec 
to proceed one time step on 1 node (32cpu) of Hitachi SR16000. Our largest calculations 
in the current report is the case of the diffusion of the 3D Gaussian packet in §4.1.11 It 
costs the memory of ~900 GB on NEC SX9 (8cpu) for the case of a fine mesh. Note that 
this specification of the computational speed is obtained within the basic optimization and 
the automatic parallelization. A parallel version of the numerical code for massive parallel 
architectures is ready for tests. In future, we would need the full coverage of the sphere 
with a high resolution by Nj. = 400, Nq = 64 and A^^ = 128, for example, for spatial grid. 
The memory requirement would be ~6 TB for the program and ~26GB for the neutrino 
distribution, which are available on the recent supercomputers. We plan to perform such 
large scale simulations after we optimize the numerical code with the parallelization. 
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z 




Fig. 1. — Definition of tlie neutrino direction. Tlie direction of neutrino propagation is 
specified by tlie neutrino angle variables, 9^ and in the spherical coordinate system. The 
neutrino angle, dy, is measured from the radial direction. The unit vectors along the radial 
and theta directions are depicted in the spherical coordinate (top). The unit vector of the 
phi direction is defined in the right-handed system (bottom). The neutrino angle, (f)^, is 
measured from the theta direction. 
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Fig. 2. — Diffusion of the Gaussian packet under spherical symmetry. The radial profiles of 
the neutrino density at the initial and final steps are shown in the top panel. The radial 
profile of the number flux at the final time step is shown in the lower panel. The solid lines 
and the cross symbols denote the analytic solution and the numerical results, respectively. 



j^Q I I I I I I I I l_J I I I I '-I I I l_J 

10 100 1000 

number of grid points 

Fig. 3. — Relative deviations of the numerical results by the 3D code for the diffusion of 
Gaussian packet under spherical symmetry as a function of the number of radial grid points. 
The relative deviations of the neutrino densities from the analytic solutions are averaged 
over all radial grid points. The dotted line shows reference proportional to iV~^. 
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Fig. 4. — Early and final distributions of neutrinos in the diffusion of 2D Gaussian packet by 
contour plot on the 2D plane. The number density in fm~^ is expressed by the color code. 
The arrows express the neutrino fluxes. 
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Fig. 5. — Neutrino distributions in the diffusion of 2D Gaussian packet. The cross symbols 
show the initial and final profiles in the numerical results along Z = —8.1 x 10^ cm (near 
the center of square). The two solid lines show the analytic solutions. 



10 100 1000 

number of grid points 



Fig. 6. — Relative deviations of the numerical results for the diffusion of 2D Gaussian packet 
as a function of the number of radial grid points. The relative deviations of neutrino densities 
from the analytic solutions are averaged over grid points along Z ~ (near the center of 
square). The dotted line shows reference proportional to N~'^. 
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Fig. 7. — Initial and final distributions of neutrinos in the diffusion of 3D Gaussian packet. 
The neutrino densities in fm~^ are shown by color contour plot on the 2D plane along 
X = 1.0 X 10^ cm. The isosurfaces of the neutrino densities (1.0 x 10"^, 1.0 x 10^^, 4.5 x 10~^ 
fm~^) are also shown by the color code. The arrows express the neutrino fluxes. 



- 69 - 




Fig. 8. — Magnitude of relative deviations in the diffusion of 3D Gaussian packet. The 
relative deviations of the neutrino densities by the numerical calculation from the analytic 
solutions are shown by contour plot on the plane along y = —1.6 x 10^ cm. 
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Fig. 9. — ID advection of a step-like distribution. The neutrino density as a function of 
radial coordinate is shown by solid lines with cross symbols for the initial time (top), the 
time at 1.0 x 10~^ s (middle), 2.0 x 10~^ s (bottom). The analytic position of the wave front 
is shown by dashed lines. 



- 71 - 



I 



6x10 

5 
4 
3 
2 
1 




6x10 



1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 


(a) : 




\ \\ 












\ \ 






\ \ 






\ '\ 






\\ 






\\ 








V 






\ \ 






\ \ J 








V ^ 



^ 



0.996 



0.998 



0.996 0.998 1.000 1.002 1.004x10** 
R [cm] 




1.000 1.002 ioo4^ioS 
R [cm] 



Fig. 10. — Neutrino densities as a function of radial coordinate for the time at 2.0 x 10"^ 
s by the calculations with different time steps and resolutions, (a) The upper panel shows 
the different cases of At=10~^ s (long-dashed line), 10"'' s (dotted line) and 10"® s (solid 
line) for the radial grid with Nr=100. (b) The lower panel shows the cases for the different 
resolutions with Nr=100, At=10"'^ s (dotted line) and N^=1000, At=10"® s (sohd line). The 
analytic position of the wave front is shown by dashed lines. 



Fig. 11. — Searchlight beam test in the 3D box. The neutrino densities are shown by color 
code in the lower panel. The source of neutrinos to one direction is set at the right boundary 
of the box. The beam propagates along the direction with a gradual spread due to the 
numerical diffusion. The numerical grid of the 3D box is shown in the upper panel. 
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Fig. 12. — (a) Energy spectra at the radial position of 98.4 km. The sohd hne and the sohd 
hne with symbol denote the spectra obtained by the formal solution and the computation 
{Nq^ = 24), respectively, (b) Relative errors of the spectra by the computation with respect 
to the formal solutions. The dotted, dot-dashed and solid lines denote the numerical results 
for Ng^ = 6, 12 and 24, respectively. 
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Fig. 13. — (a) Time evolution of the neutrino populations at the two energy grid points 
toward the equilibrium value for the dense matter in the central core, (b) Evolution of the 
neutrino spectra at selected time steps toward the equilibrium. The lines correspond to the 
time steps, 2.0xl0-^ 5.4xl0-^ 1.6x10^^ 5.2x10"^ 1.3x10^^ 3.3x10"^ sec from bottom 
to top. The solid line expresses the Fermi-Dirac distribution. 
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Fig. 14. — Radial profiles of the density, temperature and neutrino chemical potential in the 
supernova core (a) when the central density is 10^^ g/cm^ during the collap se and (b) at 100 
ms after the bounce in the spherical simulation by lSumiyoshi et al.l (l2005l ). 
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Fig. 15. — (a) Radial distributions of neutrinos for the profile at the central density of 10^^ 
g/cm^ during the collapse. The neutrino densities obtained by the 3D code are shown by 
cross symbols for three species. The densities from the spherical calculation (see the main 
text) are plotted by the solid lines, (b) Relative errors of the densities by the 3D code with 
respect to those by the ID code. 
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Fig. 16. — (a) Neutrino number fluxes for tlie profile at the central density of 10^^ g/cm^ 
during the collapse. The neutrino fluxes obtained by the 3D code are shown by cross symbols 
for three species. The neutrino fluxes from the spherical calculation (see the main text) are 
plotted by the solid lines, (b) Relative errors of the fluxes by the 3D code with respect to 
those by the ID code. 
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Fig. 17. — Same as Fig. [151 but for the profile at 100 ms after the bounce. 
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Fig. 18. — Same as Fig. [161 but for the profile at 100 ms after the bounce. 
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Fig. 19. — Effective mean free paths for tlie three species of neutrinos with the energy of 34.0 
MeV for the profile at 100 ms after the bounce. The results by the 3D code are shown by 
symbols as functions of radius. The mean free paths by the spherical calculation are shown 
by the lines. The name of neutrino reactions are indicated by the notation in §3.21 
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Fig. 20. — Ratios of the effective mean free path calculated in the 3D code with respect to 
that used in the spherical calculation for the profile at 100 ms after the bounce. The ratios 
are shown for the pair process (pap) and the nucleon-nucleon bremsstrahlung (nbr) with the 
neutrino energy of 34.0 MeV. The solid, dashed and dot-dashed lines show the ratios for Ve, 
Ve and z/^, respectively. 
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Fig. 21. — Emissivities for the three species of neutrinos for the profile at 100 ms after the 
bounce. The resuhs by the 3D code are shown by symbols as functions of radius. The 
emissivities obtained by the spherical calculation are shown by the lines. The name of 
neutrino reactions are indicated by the notation in §3.21 
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Fig. 22. — Eddington factors and flux factors as functions of radius for the profile at 100 ms 
after the bounce. The solid, dashed and dot-dashed thick lines display the results for z/g, 
and v^^ respectively. The solid, dashed and dot-dashed thin lines show the corresponding 
quantities in the spherical calculation. 
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Fig. 23. — Deviation of the Eddington factors and flux factors by the 3D code from those 
by the ID code as functions of radius for the profile at 100 ms after the bounce. The solid, 
dashed and dot-dashed lines denote quantities for z/e, Ve and z/^, respectively. 
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Fig. 24. — Luminosities and average energies as functions of radius for the profile at 100 ms 
after the bounce. The sohd, dashed and dot-dashed thick hues denote quantities for i/g, t'e 
and z/^, respectively. The corresponding quantities from the spherical calculation are shown 
by thin lines for comparison. 
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Fig. 25. — Radial profiles of the density and temperature of the deformed supernova core for 
the two directions with the polar angles oi fj, — 8.2 x 10~^ (near the equator: dashed hne) 
and fjL — 0.98 (near the pole: solid hne). 
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Fig. 26. — Contour plot of the density of the deformed supernova core with the paths of 
neutrino propagation to evaluate the formal solution at the locations denoted by symbols. 
The square and triangle symbols denote the two locations along the polar directions with 
fi = 8.2 X 10~^ and fi = 0.98, respectively. The solid lines show the paths along the 
neutrino angle, fi^ = 0.99519 for the case of Ng^ = 24. The paths along the neutrino angle. 
Hi, = 0.93247 and 0.98156, for the case of Ng^ = 6 and 12 are also shown by long-dashed 
and dashed lines for the two locations, respectively. Note that interpolation is made to plot 
smoothly contours. 
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Fig. 27. — Upper panels: energy spectra at the radial grid points 98.4 km (left) and 198 
km (right) along the polar direction /i = 8.2 x 10^^. The energy spectra evaluated by the 
3D code (solid lines with cross symbols) with Nq^ = 24 and N^^ = 6 are shown with the 
formal solutions (solid line). Lower panels: relative errors of the spectra by the 3D code 
using {Ne^, N^^) = {6,6), (12,6), (24,6) and (12,12) with respect to the formal solutions are 
shown by dotted, dot-dashed, solid and dashed lines, respectively. 
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Fig. 29. — Contour plots of the profiles in the deformed (axially symmetric) supernova core 
on the meridian slice with a constant 0. The densities of electron-type neutrinos, electron- 
type anti-neutrinos and mu-type neutrinos are shown by color codes in left top, right top 
and left bottom panels, respectively. The neutrino densities are shown here in the unit of 
fm~^=10~^^ cm~^ and in log-scale. The contour plot of the temperature is also shown in 
right bottom panel. 
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Fig. 30. — Radial profiles of the densities and fluxes for the three neutrino species in the 
deformed supernova core along the two polar directions on the constant 0-slice. The densities 
and fluxes are shown as a function of the radius in the two directions with (i — 8.2 x 10~^ 
(near the equator: dashed line) and ii = 0.98 (near the pole: solid line). 
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Fig. 31. — Contour plots of the neutrino fluxes for electron-type anti- neutrinos in the de- 
formed (axially symmetric) supernova core on the constant 0-slice. The radial and polar 
components of flux are shown by color codes in top and bottom panels, respectively. The 
fluxes are shown in the unit of fm~^s^^. 
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Fig. 32. — Radial profiles of the flux factors (top) and average energies (bottom) in the 
deformed supernova core along the two polar directions with fi = 8.2 x 10~^ (thick lines) 
and fi = 0.98 (thin lines). The solid, dashed and dot-dashed lines denote quantities for Ue, 
Vf, and z/^, respectively. 
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Fig. 33. — Elements of the Eddington tensors in the deformed supernova core as a function 
of radius along the two polar directions with fi = 8.2 x 10~^ (near the equator: top) and 
/i = 0.98 (near the pole: bottom). The diagonal elements {k^^, k^^, k"^"^) are shown by 
solid, dashed and long-dashed lines, respectively. The non-diagonal elements (A;^^, fc'^'^, k^"^) 
are shown by dot-dashed, dot-dot-dashed and dotted lines, respectively. The thin, normal 
and thick lines denote the elements for the neutrino energies, 4.89, 12.9 and 34.0 MeV, 
respectively. 
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Fig. 34. — Energy spectra by the formal (black) and steady state (red with cross symbols) 
solutions at the four locations with the radius of 98.4 km in the 3D deformed supernova core. 
Left and right panels show the spectra on the meridian slices near the zx- and yz- planes, 
respectively. The top and bottom panels correspond to the directions near the equator and 
pole, respectively. See the main text for the details of the locations. 
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Fig. 35. — Contour plots of the neutrino distributions on the meridian shce with (j) = 0.436 
radian (near zx-plane) for electron-type anti-neutrinos in the 3D deformed supernova core. 
The density, radial, polar and azimuthal components of the flux are displayed by color codes 
in left top, right top, left bottom and right bottom panels, respectively. The densities are 
shown in the unit of fm^^=10^^^ cm~^ and in log-scale. The fluxes are shown in the unit of 
fm"^s~^. 
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Fig. 36. — Same as Fig. |35l but for the meridian slice with = 1.309 radian (near yz-plane). 
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Fig. 37. — Contour plots of the angle moments, (yU^), (top) and the average energies, (e), 
(bottom) for electron-type anti-neutrinos for the slices with = 0.436 (left) and = 1.309 
(right) radian. Note that the covered scales are different between top and bottom panels. 
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Fig. 38. — Radial profiles of the elements of the Eddington tensor along the two polar 
directions with /i = 8.2 x 10^^ (top) and fi = 0.98 (bottom) on the slice with (p = 1.309 
radian (near yz-plane). The elements are shown for electron- type anti- neutrinos with the 
neutrino energy of 34.0 MeV. The diagonal elements (thick) are shown by solid, dashed and 
dot-dashed lines for k"^^, k^^, k'^'^, respectively. The non-diagonal elements (thin) are shown 
by solid, dashed and dot-dashed lines for k"^^, k^'^, k^"^, respectively. 



